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The Determination of Round-Trip Planetary 
Reconnaissance Trajectories 


RICHARD H. BATTIN* 
Massachusetts Institute of Technology 


Summary 


Some interesting geometrical and analytical properties of 
free fall trajectories are developed and subsequently exploited 
to provide a technique for determining interplanetary transfer 
orbits. After a general discussion, the determination of pro- 
pulsion-free round trips to the planet Mars is considered. Asa 
preliminary step, the reconnaissance problem is analyzed for a 
simple two-dimensional model of the solar system, assuming 
circular planetary orbits. 

Finally, a method is described for determining round-trip, non- 
stop reconnaissance trajectories in a three-dimensional model 
with elliptical planetary orbits. The results from the simplified 
approach are compared with those obtained from the true model. 
It is found that several important features of the trajectory prob- 
lem are basically three-dimensional in nature and that the simpli- 
fied model is inadequate for their description. 


(1) Geometrical and Analytical Properties 
of Trajectories 


(1.1) General Conic Trajectories 


Ww" A BOopy is in motion under the action of an 
attractive central force that varies as the inverse 
square of the distance, the path described will be a 
conic whose focus is at the center of attraction. The 
particular conic (ellipse, hyperbola, or parabola) is 
determined solely by the velocity and the distance 
from the center of force. In many problems of inter- 
planetary flight it is frequently convenient to analyze 
the flight of a space ship by considering its motion to be 
influenced by only one celestial body during any one 
period of time. Therefore, we shall begin our study of 
interplanetary trajectories by considering the purely 
geometrical problem of determining the various conic 
paths that connect two fixed points and that have a 
focus coinciding with a fixed center of force. 


Presented at the Space Flight Mechanics Session, Twenty- 
Seventh Annual Meeting, IAS, New York, January 26-29, 1959. 

* Assistant Director, Instrumentation Laboratory, Depart- 
ment of Aeronautics and Astronautics. 


There are many equivalent definitions of conics; 
however, we shall find the following ones most con- 
venient for our purposes: 


Ellipse: The locus of points the sum of whose dis- 
tances from two fixed points (foci) is constant. 
Hyperbola: The locus of points the difference of 
whose distances from two fixed points (foci) is con- 
stant. 

Parabola: The locus of points equally distant from 
a fixed point (focus) and a fixed straight line (direc- 
trix). 


The familiar elements of these conics are shown in 
Figs. 1-1, 1-2, and 1-3 and will serve to orient the reader 
during the remainder of the discussion. 

Consider now two fixed points P and Q and a center 
of force fixed at a point F. For parabolic paths these 
three points are sufficient to determine the two parabo- 
las connecting P and Q with a focus at F. However, 
for elliptic and hyperbolic paths, further specifications 
are required to determine a unique trajectory. Let us 
examine first the characteristics of the family of elliptical 
paths. 


(1.2) Elliptical Trajectories 


Consider the three points F, P, and Q in Fig. 1-4 and 
let it be required to find an ellipse with a focus at F 
that connects the two points P and Q. If the location 
of the second focus F* (sometimes called the vacant 
focus) is fixed, the problem is solved and the path is 
determined. Since the point F* cannot be placed 
arbitrarily, it will be of interest to find the locus ot the 
foci of all ellipses that satisfy the conditions of the 
problem. 

To this end, let us designate the line segments in 
Fig. 1-4 as follows: 


FP=n, FQ=n, PQ=c (1.2.1) 
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MAJOR AXIS ~ 20 
MINOR AXIS ~ 2b 
LATUS RECTUM - 2 
ECCENTRICITY € <1 
PF . PF” 20 


V1-(b 0)? 
Le of] 


MAJOR AXIS = 20 
MINOR AXIS = 2b 

LATUS RECTUM - 2. 
ECCENTRICITY -€>1 

SLOPES OF ASYMPTOTES ~ bo 
PF” -PF 20 


P*F -P°F* ~20 
€= V1+(b/0)? 
t= ole? 1) 


DISTANCE FROM VERTEX TO FOCUS - 0 
LATUS RECTUM ~ 2£ 


PF - PN 


20 = 
Fic. 1-1 (Top). Ellipse. Fic. 1-2 (Center). Hyperbola. 
Fic. 1-3 (Bottom). Parabola. 


(In our discussion we shall assume 72 > 7; obvious 
changes in the results can be made if the reverse in- 
equality holds. The case for which 7; = 12 is quite 
special and, indeed, almost trivial.) Since P and Q 
must both lie on the ellipse, the point F* must be se- 
lected such that 


PF* + PF = QF* + OF = 2a 
or, equivalently, 
PF* =2a-—1n, QF* 


Thus, for an ellipse of major axis 2a, the point F* is de- 
termined as the point of intersection of two circles 
centered at P and Q with respective radii 2a — 7, and 
2a — 7. A number of such circles have been con- 
structed in Fig. 1-4 for different values of the major 
axis 2a. 

We may make at once several interesting obser- 
vations: 

(a) If the selected value of 2a is too small, the circles 
will not intersect. Thus, there is a smallest value, 
2am, below which no elliptical path is possible. When 
a@ = Gm, the two circles are tangent and the point of 
tangency F,,* lies on the line PQ. Thus, a, may be 
determined from 

(2am — + (2am — n) = 
to give a value 
2am = (1) + re + c)/2 (1.2.2) 


which is just one-half of the perimeter of the triangle 
FPQ. For later convenience, we introduce the nota- 
tion 


(1.2.3) 


so that 
2am = 
The point F,,* divides the line PQ in such a way that 
PF,* =s—1n, QF,»* =s— fe 


(b) When a > a,, the pair of circles intersect in two 
points F* and F*. Thus, there are, in general, two 
different elliptical paths connecting P and Q having 
the same length major axis but with vacant foci equi- 
distant from and falling on opposite sides of the line 
PQ. For any value of 2a, the focus F* is a greater dis- 
tance from F than the corresponding conjugate focus 
F*, Therefore, the ellipse with focus at F* hasa larger 
eccentricity and a smaller latus rectum than the ellipse 
with the same length major axis and focus at F*. 

As a point of interest, the values of a; and a» used in 
the construction of Fig. 1-4 are 


2a, = 2am + r,/2, 2a2 = 24m + 


(c) Each vacant focus is so located that the distances 
from P and Q have a constant difference 7 — 1. 
Thus, the locus of these foci is a hyperbola; P and Q are 
the foci of the hyperbola; re — 7 is the length of its 
major axis; and c/(r2 — 7) is the eccentricity. The 
asymptotes of the hyperbola have slopes given by 


Slopes of the asymptotes = 


+ [2V — sin (0/2) (1.2.4) 


The second form of Eq. (1.2.4) shows explicitly how the 
hyperbolic locus of F* varies with the angle between 
FP and FQ. 
Minimum-Energy Ellipse 

The point F,,* defines the so-called ‘‘minimum- 
energy” elliptical path from P to Q. The kinetic 
energy of a body, moving in free fall along an elliptical 
arc, is proportional to [(1/r) — (1/2a)], where r is the 
distance from the center of force. Thus, the kinetic 
energy at any point is a minimum when the major axis 
of the path has the minimum value of 2a,,. We may 
determine the semi-latus rectum /,, and the eccentricity 
€, of the minimum energy path in the following manner: 


Since = 
we have, from Fig. 1-4, 


(2€mOm)? = [(s — rm) sin Z PQF]? + 
[r2 — (s — m) cos Z POF}? 


Using the trigonometric identity 
cos Z POF = [2s(s — n)/rec] — 1 
we have 


= s? — 4s(s — 1) (Ss — 


But, for an ellipse, in general, 


(2ea)? = 4a(a — /) 
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so that (ZenBm)? = s? — Qsly 


Hence, 


= — n)(s — 72)/c = 
(mr2/c)(1 — cos 6) (1.2.5) 


The eccentricity €, may then be found from 


In = (s Em”) 


Osculating Ellipse 

The point Q coincides with the point of aphelion for 
the ellipse whose focus F,* is determined as the inter- 
section of the hyperbolic locus and the line QF. All 
elliptical paths from P to Q with foci to the right of the 
line QF reach aphelion before arriving at Q, while those 
with foci on the left reach aphelion after passing 
through Q. Since the path is tangent at Q to a circle 
of radius 7 and centered at F, we shall use the term 
“osculating ellipse’ when referring to this trajectory. 
We may again employ simple geometrical arguments 
to determine the major axis 2d), eccentricity «, and 
latus rectum Jy. 

Referring to Fig. 1-4, we have, from the definition 
of an ellipse, 


QF,* + QF = PF + PF,* = 2a 
However, since Q, F, F,* are colinear, then 
QFy* = re — edo = 2a — fr 
Hence, PF,* = 2re — — 11 


Considering now the triangle FPFy)*, we use a funda- 
mental trigonometric identity to obtain 


re[r2 — (2re — — 1 + cos 6 
which may be solved for 2¢a@) to yield 
= 27re(r2 — 1)/[2re — r1(1 + cos 6) | 
But 2ay = — 
so that 


Qe = {1 + [(r2 — n)/(re — n cos (1.2.6) 


The corresponding value of the semi-latus rectum is 
readily found to be 


= — cos (1 — cos (1.2.7) 


Finally, by comparing Eqs. (1.2.5) and (1.2.7), we note 
the following interesting relation: 


ln = cos Z PQF (1.2.8) 


Symmetrical Ellipse 

To complete our discussion of elliptical trajectories, 
let us consider the ellipse whose focus F,* is determined 
as the intersection of the hyperbolic locus and the line 
through F that is parallel to PQ. This path is of some 
interest because of its symmetry—i.e., the point P 
bears the same relationship to perihelion as the point 
Q does to aphelion. The elements of this ellipse, 2a,, 
e, and /,, are readily obtained since the polygon 
PQF,*F is an isosceles trapezoid. 


We find 2a4,=nHn+Pre (1.2.9) 


and 


= [nre(r1 + r2)/c?] (1 — cos @) = 
+ (1.2.10) 


(1.3) Hyperbolic Trajectories 


Consider again the three points F, P, and Q shown 
now in Fig. 1-5 and let it be required to find a hyper- 
bola with focus at F that connects the two points P 
and Q. Again the problem is solved if the location of 
the second focus F* is determined. However, since the 
point F is an attractive focus, we are interested only 
in the branch of the hyperbolic path that is concave 
with respect to F. 

The points P and Q must both lie on the concave 
branch of the hyperbola, so that the point F* must be 
selected such that 


PF* — PF = QF* — QF = 2a 


Thus, for a hyperbola whose major axis is 2a, the point 
F* is determined as the point of intersection of two 
circles centered at P and Q, with respective radii 
2a + m, and 2a + ro. Three such pairs of circles have 


Fic. 1-4. Locus of the vacant foci for elliptical trajectories. 


Fic. 1-5, Locus of the vacant foci for hyperbolic trajectories. 
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Fic. 1-6. Parabolic trajectories. 


been constructed in Fig. 1-5 for values of a = 0, a, 
and ds, where a; and ad» where chosen as 


2a; = 7 = 7) 


We make the following interesting observations: 

(a) All points of intersection of the circle pairs fall 
outside of the circle centered at P and of radius7;. One 
may readily verify that, for hyperbolic paths from P 
to Q that are convex with respect to the focus F, the 
vacant foci F* all lie within this circle. 

(b) The pair of circles intersect in two points F* 
and F*, so that there are two different hyperbolic paths 
connecting P and Q having the same length major axis 
and with vacant foci equidistant from and falling on 
opposite sides of the line PQ. For any value of 2a, the 
hyperbola with focus at F* has a larger eccentricity 
and a larger latus rectum than the corresponding 
hyperbola with focus at F.* 

(c) Each vacant focus is so located that the difference 
of its distances from Q and P is r, — r,.. Thus, the 
locus of these foci is the conjugate branch of the hyper- 
bolic locus of the foci for the elliptical paths considered 
previously. 

The foci F)* and /,*, corresponding to a zero-length 
major axis, are extreme cases in that infinite velocities 
are required to describe the associated paths. The 
path with vacant focus at F;* is the straight line from 
P to Q;i.e., the hyperbola for which a = Oande = ~. 
Corresponding to the focus, /)*, the path is composed 
of the two straight-line segments from P to F and from 
F to Q and is the hyperbola for which a = 0 and e€ = 
sec (6/2). 


(1.4) Parabolic Trajectories 


There are two parabolic paths with focus at F that 
connect the two points P and Q. To determine the 
axes of these parabolas, we shall first locate their direc- 
trices. Referring to Fig. 1-6 and recalling the defini- 
tion of a parabola, the directrices D,D;’ and DD,’ 
are obtained as the two common tangents of the two 
circles centered at P and Q with respective radii 7; and 
ro. The axes AA,’ and A,A,’ of the two parabolas are 


the normals to the corresponding directrices through 
the focus F. The vertices 1; and I, are the mid- 
points of that portion of the axes included between the 
focus F and the directrices. By elementary geometry 
one can show that the axes A,A,’ and AA,’ are parallel 
to the asymptotes of the hyperbolic locus of the vacant 
toci of the elliptical and hyperbolic paths considered 
previously. Thus, Eq. (1.2.4) gives the slopes of these 
axes with respect to the line PQ. 

Again, it is an elementary exercise in geometry to 
determine for each parabola the semi-latus rectum or, 
equivalently, the distance from focus to directrix. One 
obtains 


= 2FV; = [4(s — — X 
[Vs/2+ V(s — 0/2)? (14.1) 
l, = = [4(s — re) X 


[Vs/2 — V(s — 0)/2}2 (1.4.2) 


(1.5) Latus Rectum and Eccentricity 


Although we have found it possible to obtain the 
jatus rectum and eccentricity of several special conic 
paths by employing geometrical arguments, it would be 
impractical to use similar techniques for the general 
case. Instead, let us tackle the problem of obtaining 
an analytic functional relationship between the latus 
rectum and major axis of the conic by using the polar- 
coordinate equation 


r = 1/(1 + € cos ¢) (1.5.1) 


where ry and @¢ are polar coordinates of a point on the 
conic whose focus is at the origin. The conic is either 
an ellipse or a hyperbola, according to whether the 
eccentricity ¢ is less than or greater than one, respec- 
tively. The semi-latus rectum / is related to the eccen- 
tricity « and the semi-major axis a by 


/=a(1 — &) (ellipse) 
l = a(ée? — 1) 


Since the points P and Q of Fig. 1-7 both lie on the 
conic, then, from Eq. (1.5.1), we have 


(1.5.2) 


(hyperbola) (1.5.3) 


Fic. 1-7. General conic trajectory. 
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cos (y + 6) = — 1 
(1.5.4) 


Consider the trigonometric identity 


cos? (y + 6) — 2 cos (y + 8) cos y cos 6 + 
cos? y — sin? @ = 0 


and substitute from Egs. (1.5.1)--(1.5.3) to obtain 


— re)? — Zaryre(l — re) (1 — cos 6 + 
— — (a ¥ sin? @ = 0 


where the selection of upper or lower signs in the last 
term is made, respectively, according to whether the 
conic is an ellipse or a hyperbola. If we now collect 
terms in powers of / (note that the coefficient of /? is 
simply ac?) and make further obvious simplifications, 
we find 


acl? — ryre(1 — cos 6) [2a(r; + re) 


nre(1 + cos @)|/ + ar,*r.2(1 — cos 0)? = 0 (1.5.5) 


Introducing the symbol s, which is defined by Eq. 
(1.2.3), permits the bracketed expression of Eq. (1.5.5) 
to be rewritten in the form 


2a(r, + re) F + cos = 2a(2s — c) 


2s(s — c) = ¥2s(s — c F 2a) — 2ac (1.5.6) 


In the remaining part of the derivation we shall con- 
sider separately the ellipse and the hyperbola. 
For the elliptical case, it is convenient to introduce 
the notation 
= / 
sin(a/2) = + = Vs/2a (1.5.7) 
sin (8/2) = V(r, rs 4a = Vis — ¢)/2a 
(1.5.8) 
so that we may write 


s — ¢ — 2a = —2a cos? (8/2) 
s = 2a sin? (a/2) 
c = 2a [sin® (a/2) — sin? (8/2)] 


The right-hand side of Eq. (1.5.6) may then be written 
as 


Sa? sin? (a@/2) cos? (8/2) — 4a*[sin? (a/2) — 
sin? (8/2)] = 4a? sin? (a/2) cos 8 + 4a? sin® (8/2) = 
2a*{ sin? [(a + B)/2] + sin? [(a — 8)/2}} 
Substituting this into Eq. (1.5.5) and further noting 
that 
— cos 0) = 2(s — n)(s — 
gives 
acl? — 4a2(s — n)(s — sin? [(a + B)/2] + 
sin? [(a — B)/2]} 1 + 4a(s — n)°%(s — = 0 
Finally, multiplying through by c?/a and using the 
identity 
c = 2asin [(a + B)/2] sin [(a@ — 8)/2] 


yields 


cil? — 4a(s — r)(s — re)} sin? [((a + B)/2] + 
sin? [(a — )/2]} cl + 16a%(s — (s — X 
sin? [(a@ + 8)/2] sin? [(a — 8)/2] = 0 


from which the two roots 


l = [4a(s — — re)/c?] X 


sin? [((a + 8)/2] (ellipse) (1.5.9) 
are obtained. 
For the hyperbolic case, an analogous procedure is 


followed after introducing the notation 


sinh (a/2) = W(n + r+ 0)/4a = Vs (1.5.10) 


Vin + c) 4g = V(s c)/2a 


sinh (8/2) 


We may then write 
s — c + 2a = 2a cosh? (8/2) 
s = 2a sinh? (a/2) 
c = 2a [sinh? (a/2) — sinh? (6/2) | 
and the right-hand side of Eq. (1.5.6) becomes 
2a? { sinh? [(a + 8)/2] + sinh? [(a@ — 8)/2}} 
Substituting into Eq. (1.5.5), using the identity 
c = 2a sinh [(a@ + 8)/2] sinh [(a@ — 8) 2] 


and following the same steps as for the elliptical case, we 
obtain 


= [4a(s — n)(s — r:)/c?] X 


sin'i? + B)/2] (hyperbola) (1.5.12) 


Our analytical and geometrical results are consistent 

that is, to each value of a, there correspond two el- 
lipses or two hyperbolas that satisfy the requirements 
of the problem. Furthermore, in the elliptical case 
the conditions for minimum energy, as expressed by 
Eq. (1.2.2), imply a = 7 in Eq. (1.5.7). Also, since 

it follows that 
sin? [(@ + 8)/2] > sin® [(a — 8)/2] 

Therefore, if /, and /_ are used to denote the two 
roots given by Eq. (1.5.9), in accordance with the 
particular choice of sign, we have 


Since the distance L between the foci of an ellipse may 
be written as 
L = 2Va(a — J) 
L, L. 


with the subscripts having the obvious interpretation. 
Therefore, in Eq. (1.5.9), the upper sign corresponds 
to the focus F*, shown in Fig. 1-4, while the lower sign 
corresponds to the focus F*. 

Similarly, in the hyperbolic case, 


0<B<a 


it follows that 


5.1) 
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5.2) 
5.3) 
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so that 

sinh? [(a + 8)/2]| > sinh? [(a — 8)/2] 

Again denoting by /, and /_ the two roots given by Eq. 
(1.5.12), we have 


as before. On the other hand, the distance L between 
the foci of a hyperbola is 


L = 2Va(a + 


so that iy 


Therefore, in Eq. (1.5.12), the upper sign corresponds 
to the focus F*, shown in Fig. 1-5, while the lower sign 
corresponds to the focus F*. 

We may check our general formula for the latus 
rectum with the special results obtained previously by 
geometrical arguments. For example, in the mini- 
mum-energy case, where a,, = s/2, we have 


2, V(s — o/s 


Substitution into Eq. (1.5.9) produces the same ex- 
pression for /,, as obtained in Eq. (1.2.5). The two roots 
normally obtained are identical, showing that only one 
minimum-energy path exists from P to Q. 

For the case of the symmetrical ellipse in which, from 
Eq. (1.2.9), 


a, = (mn + 12)/2 = (2s — c)/2 


we see at once that the corresponding values of a and 8 
are such that 


a,+p,=7 


Eq. (1.2.10) for /, follows immediately from the general 
expression for / if we choose the upper sign in Eq. (1.5.9). 
We shall close this discussion by making one final 
observation. By using Eqs. (1.5.7) and (1.5.8), we 
may expand Eq. (1.5.9) as 


l= [4(s re) x 
{Vs/2V1 — [(s — c)/2a| + 


V1 — (s/2a) V(s — 0/2}? 


and determine the limit of / as a becomes infinite. Re- 
ferring to Eqs. (1.4.1) and (1.4.2), we see at once that 
lim/, =), lim/_ 

a—> oa 

which shows that the parabolic paths from P to Q are 
the limiting forms of elliptical paths with infinite major 
axes. A similar expansion of Eq. (1.5.12) and subse- 
quent determination of the limit for increasing a also 
shows these same parabolas as the limiting forms of 
hyperbolas with infinite major axes. 

Finally, one may consider, from an analytic point 
of view, the shape of the limiting hyperbolas as a ap- 
proaches zero. By calculating the limiting slope of 
the asymptotes 


lim (6?/a?) = lim (//a) 
a—0 a—0 
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the results stated at the end of Section (1.3) may be 
verified. 


(1.6) Time of Flight for a Conic Trajectory 


There are many interesting and sometimes surprising 
properties of conic trajectories. For example, one 
would scarcely anticipate that the period P of elliptic 
motion would depend only upon the major axis of the 
ellipse and not at all upon the eccentricity. The pre- 
cise relationship is 


P = Vai/u (1.6.1) 


where a is the semi-major axis of the ellipse and uy is a 
constant of such a value that u/r? gives the magnitude 
of the force at a distance r from the center of force. 
The fact that a particle moving in free fall along a 
conic trajectory has a velocity magnitude |’ at any 
point that is a function only of the distance r from the 
center of force and the major axis is another example 
of one of the more interesting properties of conics. 
Here, again, nc dependence upon eccentricity is in- 
volved and the velocity is determined accordingly by 


V? = p[(2/r) — (1/a)| (ellipse) (i.6.2) 
V? = p[(2/r) + (1/a)] (hyperbola) (1.6.3) 
V2 = 2u/r (parabola) (1.6.4) 


Perhaps the most remarkable theorem in this connec- 
tion was the one originally discovered by Lambert, and 
subsequently proved analytically by Lagrange, having 
to do with the time to traverse an elliptic are under 
conditions of free fall. Lambert showed that this time 
depends only upon the length of the major axis, the 
sum of the distances of the initial and final points of 
the arc from the center of force, and the length of the 
chord joining these points. (Actually, the theorem is 
true for a general conic.) In terms of our notation, 
if T is the time to describe the arc from P to Q shown 
in Fig. 1-7, then Lambert’s theorem states that 


T = T(a, nm + fe, c) (1.6.5) 


We are again astonished to find that the ellipticity is 
not involved. 

In this paper we shall not be concerned with a demon- 
stration of the truth of Lambert’s result, but will in- 
stead exploit an idea suggested by Lagrange to arrive 
at the precise analytical form of the functional relation- 
ship implied in Eq. (1.6.5). Since we have already seen 
that there are two paths from P to Q for each conic 
having the same values of a, 7, 72, and c, we must expect 
to find two different analytical expressions—one valid 
for each of the two paths.t 

Consider first an elliptical arc from P to Q and, more 


+ For an ellipse connecting the points P and Q, there is, of 
course, a choice between two paths, according to whether the 
ellipse is traversed in the clockwise or counterclockwise direction. 
In our calculation of the time of flight 7, we shall always assume 
counterclockwise motions from P to Q. The time to traverse the 
clockwise path may be obtained by subtracting T from the total 
period. 
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specifically, one whose vacant focus F* lies along the 
lower branch of the hyperbolic locus shown in Fig. 1-4. 
Since the time 7 to traverse this are does not depend 
on the eccentricity, we may consider instead the arc of 
a very flat ellipse, obtained by letting « approach unity 
in such a way that a, r; + 72, and c remain unchanged. 
Then, according to Eq. (1.6.5), the time to traverse 
the original are subtended by the chord c and the time 
to traverse the corresponding arc of the flattened 
ellipse are the same. In this limiting process, as « tends 
to unity with a fixed, the foci move out to the extremi- 
ties of the ellipse and the entire curve flattens out to 
coincide with the major axis in the limit. Various 
stages, as the limit is approached, are shown in Fig. 1-8. 
The same time is required to traverse each of the three 
elliptical arcs shown in the Figure. Although the 
straight-line distance from P to Q remains constant, 
the distances of P and Q from the attractive focus will 
change, but in such a way that PF + QF is invariant. 

In the limit, the trajectory is rectilinear, the arc in 
question coincides exactly with the chord c, and we 
may compute the time 7 by elementary methods. To 
this end, we use Eq. (1.6.2) in the form 


V2 = (dr/dt)? = u[(2/r) — (1/a)] 


since now all motion takes place along a straight-line 
path. Rewriting the preceding equation, we have 


dt = (1 iV [ry dr 2r — (r?/a)| 
From Fig. 1-8(c), the end points of the are P and Q are 
so located that 
PF = (n +72 
OF = (n +r+0)/2=s 
Hence, 


T = (1/Vu) — (r2/a)] 


Introducing the auxiliary quantities a and 8, defined 
by Eqs. (1.5.7) and (1.5.8), and making the following 
change in the variable of integration: 


r = a(l — cos 


we have, finally, 


T = (1 — cos 
B 
which becomes 


T = (P/2x)[(a — sin a) — 
(8 — sin 8)| 


after performing the simple quadrature and introducing 
the period P defined in Eq. (1.6.1). 

We now derive a corresponding expression for the 
time 7 when the elliptical arc from P to Q has its vacant 
focus F* along the upper branch of the hyperbolic locus 
in Fig. 1-4. In this case, the chord c intersects the line 


(ellipse) (1.6.6) 


connecting the two foci F and F* and this characteristic 
must be preserved when we pass to the limit of a very 
flat ellipse. This situation is illustrated in Fig. 1-9 


PF + QF 


(b) 
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Fic. 1-8. Illustration of the significance of Lambert’s theorem. 


for various stages as ¢ tends to unity. In the limit, 
the time 7 is then identical to the time to traverse the 
flattened ellipse from P to F* and from F* back to Q. 
Thus, to obtain 7, we simply add to T twice the time 
to traverse the distance F*Q; i.e., 


= 2a 
T=T+Q vo f [r dr V2r (9? a) | 


= T+ 2Vai/p (1 — cos 


Hence, 


T = P — (P 2n)[(a — sin a) + 
(8 — sin 8)} 


For the case of the minimum-energy path, Eqs. (1.6.6) 
and (1.6.7) produce identical results. We have 


(Pm = (Bm — sin Bm) | 


(ellipse) (1.6.7) 


where 


P,, = 33 ‘Qu, sin (8,,/2) = V(s c) ‘Ss 


For the symmetrical ellipse, the time-of-flight expres- 
sion is particularly simple. We have 


T, = (P,/22)(2a, — = — 26;) 


where 
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Fic, 1-9. 


P, = (1 + 
‘ cos Bs = —cosa, = c/(m + re) 


Exactly the same technique may be used to obtain 
an analytic expression for the time to traverse a hyper- 
bolic arc. Consider first a hyperbola connecting the 
points P and Q whose vacant focus F* lies along the 
upper branch of the hyperbolic locus shown in Fig. 1-5. 
To compute the time 7 to traverse this arc, we again 
let « approach unity in such a way that a, 7 + 7, and 
c remain unchanged. In the limit, when the trajectory 
is rectilinear, we may use Eq. (1.6.3) to obtain 


T= W/V u) + 


lf we introduce the auxiliary quantities a and 8, defined 
by Eqs. (1.5.10) and (1.5.11), and make the following 
change in the variable of integration: 


r = a(cosh y — 1) 
we find 


T = [(sinh a — a) — (sinh — 8)] 


(hyperbola) (1.6.8) 


The hyperbolic trajectory, whose vacant focus F* 
lies along the lower branch of the hyperbolic locus in 


Illustration of the significance of Lambe rt’s theorem. 
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Fig. 1-5, is characterized by the fact that the center 
of force F is contained in the area bounded by the chord 
cand the are that it subtends. This characteristic must 
be preserved when passing to the limit. The time 7 
is then computed as the time to traverse the flattened 
hyperbola from P to F and from F to Q; i., 


T=T+(2, V'u) dr/V 2r + a) | 


Hence, 


T = Vai uw [(sinh a — a) + (sinh B — 8B)| 


(hyperbola) (1.6.9) 


We have seen, in our earlier discussions, that the 
parabolic trajectories from P to Q may be considered 
as limiting cases of either elliptic or hyperbolic tra- 
jectories with infinite major axes. If we compute the 
limit of T and T as a becomes infinite, we should obtain 
the times 7) and 7» to traverse the parabolic arcs shown 


in Fig. 1-6. Indeed, from Eqs. (1.6.8) and (1.6.9) we 
find 
T, = lim T = (1/3) V2/u [s¥2 — (s — 07] 

(parabola) (1.6.10) 
= lim 7 = (1/3) [s*2 + (5 — 

(parabola) (1.6.11) 


The same expression for 7; in Eq. (1.6.10) may also 
be obtained using Eq. (1.6.6). On the other hand, the 
parabola PV.Q of Fig. 1-6 corresponds, in the limit of 
increasing a, to the lower branch of the elliptical arc 
with vacant focus F*. Thus, to produce the formula for 
T. from the time-of-flight expression for an elliptical 
arc, we must use the complementary form of Eq. (1.6.7) 
obtained by subtracting 7 of that equation from the 
period P. 


(2) Interplanetary Trajectories in a Simple 
Model of the Solar System 


(2.1) Departure Velocity From a Circular Orbit 


Of fundamental importance to the problem of plane- 
tary reconnaissance is the impulse in velocity needed 
in order that a spaceship may attain a suitable inter- 
planetary orbit. In the present section we shall ana- 
lyze these velocity requirements using a simplified 
model of the solar system. The basic assumptions 
will be that the planets describe circular orbits around 
the sun and that the planes of these orbits are situated 
in the plane of the ecliptic. Later, in Section (3), we 
shall consider the more complex mode] in which the 
planetary orbits are elliptical and are inclined with 
respect to the plane of the ecliptic. Furthermore, 
we shall not consider here the problem of launching 
a vehicle from the surface of a planet. We assume, 


instead, that our interplanetary voyage begins at a 
point that is sufficiently far removed from the gravita- 
tional field of the planet that we need consider only the 
The vehicle or spaceship has 


attraction of the sun. 
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an initial velocity that is the same as the orbital ve- 
locity of the planet of departure. Our problem will be 
to determine the additional velocity required to place 
the spaceship in an orbit that intersects the orbit of 
the destination planet at a predetermined point in space. 
As the first step in our analysis, we will derive ex- 
pressions for the polar-coordinate components of the 
velocity of a vehicle moving in a conic trajectory. 
From the polar equation for a general conic, 
r = 1/(1 + ecos ¢) (2.1.3) 

we obtain 

dr/dt = [esin + cos ¢)| r(d@ dl) 


from which it follows that 


(dr + [r(dd/dt) |? = 
[ee — 1 + 21 2) [r(do |? 


—(l/a) (ellipse) 
But, 0 (parabola) 
l/a (hyperbola) 


—(u a) (ellipse) 
(dr/dt)? + [r(d@/dt) |? = (2u/r) + 0 (parabola) 
(hyperbola) 


Hence, 
[r(d@/dt) |? = ul /r? (3.23) 
\~ (u/a) (ellipse) 


(dr ‘dt)? = p[(2r — J)/r?| + 0 (parabola) 
p/a (hyperbola) 
(2.1.3) 


In order to keep the energy requirements within rea- 
sonable bounds, we shall restrict our discussion of de- 
parture velocities for orbital transfer missions to those 
required for elliptical trajectories. 

Refer to Fig. 2-1 and consider a spaceship in a circular 
orbit of radius 7, about a center of attraction. The 
orbital velocity Vy is 


Vo? = 2.1.4) 


When the vehicle is at the point P, a velocity increment 
Vx is applied in such a way that the ship will leave its 
present orbit at P with velocity Vp and move along an 
elliptical trajectory to a point Q at a distance r, from 
the center of force. We wish to determine the velocity 
increment Vz as a function of the heliocentric angle 4 
through which the ship moves in its voyage from P to Q. 
In particular, we are interested in the minimum velocity 
required for any given mission. 

After the velocity increment Vp is applied, the ship 
will have a velocity Vp whose polar components 
(dr /dt)p and (rd¢@/dt)p are obtained from Eqs. (2.1.2) 
and (2.1.3). Since Vz is defined by 


Va? = (dr/dt)p? + { [r(db/dt)]» — Vo}’ 
we may use Eq. (2.1.4) to obtain 
Ve? = — (n/a) — 2Vi/n] 


(2.1.5) 
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Velocity requirements for an orbit to orbit vovage 


Let us introduce the dimensionless quantity 
AE = VR? Vo? 
which is a measure of the amount of energy needed at 
point P to transfer from a circular orbit to an elliptical 


orbit for a voyage to Q. Then, from Eq. (1.5.9) and 
a simple trigonometric identity, we have 


AE = 3 — (n/a) — 4(V ar 


sin (6/2) sin [((a + 8)/2] (2.1.6) 


Using the definitions of a and 8 given by Eqs. (1.5.7) 
and (1.5.8), we may write AF in an alternate form as 

AE = 3 — (n/a) (nV c) X 
sind{V[1/(s — c)] — (1/2a) + 


V(1/s) —(1/2a)} (2.1.7) 


From an examination of the derivative 


27, Vil (s — c)]— (1/2a) 


l 
Vi 5) | 


da a 
and Eq. (2.1.7), we note the following: 

(a) AE is a double-valued function of a having an 
infinite slope at a = ad, = s/2, the smallest value of a 
for which an elliptical path from P to Q is possible. 

(b) If we denote by AE, and AE_ the two branches 
of the function AE corresponding respectively to the 
upper and lower signs in Eq. (2.1.7), we see that 


AE, < AE_ 


with equality obtaining only when a = ad, 
(c) Asa increases, AE approaches asymptotically the 
values 


3- (rxV 2r,/c) sin 6{(1 Vs ¢) + (1 V's)] 


(d) The slope of the upper branch AE_ is always posi- 
tive, while the slope of the lower branch AE, is negative 
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Fic. 2-2. Earth to Mars departure velocity vs. semi-major axis 
of transfer ellipse for @ = 120°. 


for a near d,. The possibility that AE, will have a 
minimum for a finite value of a will be considered in the 
following section. 

A graph of the function ~/AE versus a is provided 
in Fig. 2-2. The selected mission is a voyage from 
Earth to Mars, so that 7; and rz are chosen, respectively, 
as the mean distances of the planets from the sun. 
The heliocentric angle @ traveled was taken as 120° for 
purposes of illustration. The ordinate Vpz/V,% shows 
directly the additional velocity that must be provided 
for the trip as a fraction of the Earth's orbital velocity. 
Since the Earth’s velocity is almost 100,000 ft. per 
sec., the approximate conversion of Vg into units of 
feet per second is immediate. The abscissa a is shown 
as a multiple of the astronomical unit (A.U.). 


(2.2) Minimum Departure Velocity From a Circular Orbit 


The minimum-energy trajectory discussed in Section 
(1) should not be confused with the trajectory requiring 
the minimum departure velocity. Indeed, for the 
particular case illustrated in Fig. 2-2, the additional 
velocity in excess of the orbital velocity needed to travel 
the minimum-energy path (a = d,) is about double 
the minimum departure velocity. As a matter of 
fact, the only trajectory that minimizes both the ve- 
locity relative to the sun and the velocity relative to 
the planet of departure is the so-called Hohmann ellipse. 
The Hohmann orbit is distinguished as being co-tan- 
gential to the orbit of the planet of departure and to the 
orbit of the destination planet. 

The energy requirements for the Hohmann path 
cannot be obtained directly from Eq. (2.1.7). The 
polar angle 6 is equal to m and Eq. (2.1.7) is indeter- 
minate for this value. However, because of the sim- 
plicity of the geometry, we may evaluate the elements 
of the Hohmann ellipse directly and then use Eq. 
(2.1.5) to determine the energy requirements AF. 

The points of departure and arrival for a Hohmann 
trajectory correspond, respectively, to perihelion and 
aphelion. Thus, if ay, /,, and ey, are, respectively, the 
semi-major axis, semi-latus rectum and eccentricity, 
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we have 


ay(1 — €x) =n, Qy(1 + ex) = fr 


so thatt 
= + ro) ly = 2[nre + ro) | 
Finally, from Eq. (2.1.5), we have 
AEy = 3 — [2n/(n — 2V2r/(ni +72) (2.2.1) 


The Hohmann trajectory is not necessarily the ideal 
one for interplanetary travel. The precise orientation 
of the departure and destination planets that is re- 
quired for the Hohmann path is one obvious disad- 
vantage. Another, to be discussed in Section (3), re- 
sults from the three-dimensional nature of the solar 
system. Its chief interest, therefore, lies in the fact 
that it provides, in our simple model, a lower bound for 
the energy requirements of any mission—i.e., for a 
fixed r; and r2 we have 


AEy < AE 


with equality obtaining only for @ = mw anda = (nm + 
ro) 2. 

Obviously, in planning an interplanetary voyage, it 
is of advantage to be as flexible as possible with regard 
to the time of departure. Therefore, it is of interest 
to examine the energy requirements AF for an arbi- 
trary configuration of the departure and destination 
planets. To this end, let us determine, for a fixed 7, 
72, and 6, under what conditions AE, as defined in Eq. 
(2.1.7). has a minimum value. 

It was noted in the previous section that the slope 
of the lower branch AF, of the AE curve is negative 
for values of a near ad». Referring to Eq. (2.1.8), 
we see that AL, will have a positive slope if a value of 
a exists such that 


{1/V/[1/(s — — (1/2a)} + 
(1/W(1/s) — (1/2a)] < sin 8 


Clearly, the difference between the lett-hand side of 
the preceding inequality and Vs — ¢c + Vs can be 
made arbitrarily small by choosing a to be sufficiently 
large. Thus, we are led to determine the conditions 
under which the inequality 


Vs < 27; ‘ro sin 6 
holds. For this purpose, we note that 
(Vs —c + Vs)? = 2V cos (0/2) + 
< (Vn + Vn)? 
and (2c Vv 2r, /ro sin 0)? > 8r; 


Hence, a sufficient condition for E, to possess a mini- 
mum is that 


t If we compare these results with Eqs. (1.2.2), (1.2.5), (1.2.6), 
(1.2.7), (1.2.9), and (1.2.10), we see that the Hohmann orbit, 
the minimum-energy ellipse, the osculating ellipse, and the sym- 
metrical ellipse are all identical for 6 = zx. 
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Vr < 2V2n 


If 72 < , the condition of Eq. (2.2.3) holds. There- 
fore, an elliptical path, requiring minimum departure 
velocity, to an inner planet always exists. Further- 
more, one sees that this condition holds also for a trip 
from Earth to Mars. On the other hand, for the re- 
mainder of the outer planets, the inequality does not 
hold and one might expect to find values of @ for which 
the minimum departure velocity trajectory is parabolic. 
Asa matter of fact, this rather surprising condition does 
prevail for Jupiter and the planets beyond. One 


(2.2.3) 


can always be satisfied if sin @ is made small enough. 
Thus, for the outer planets there are sectors near 6 = 0) 
and @ = a for which minimum-departure-velocity 
elliptical trajectories exist. However, the farther from 
the sun, the smaller these sectors become. 

The value of a for which AF, attains its minimum 
value (assuming that one exists) is obtained as the root 
of the equation 


Vv [1/(s _ 2a)} + 
— (1/2a)] = sin (2.2.4) 


which, after reduction to a normal form, is found to be 
of fourth degree in a. For practical purposes, the solu- 
tion of Eq. (2.2.4) is obtained more easily as follows. 
First, rewrite the equation in the equivalent form 


{1 s) 2a) + [c/s(s —c)]} + 


{1 iV (1 /s) — (1 2a) | = 21, sin 0 


and then introduce a new variable ¢ defined by 
tang = V(1/s) — (1/2a)/Ve/s(s — 0), 
(2.2.5) 
so that Eq. (2.2.4) becomes 
cos £ + cot ¢ = (2c/re sin 6) V 2¢7r; s(s — ©) 
or, alternately, 


cos + cot = 4(c/re)*/?/sin 6 Vi+ cos 8 (2.2.6) 


Eq. (2.2.6) may be solved for ¢ almost by inspection, 
using a table of trigonometric functions. With this 
value of ¢, Eq. (2.2.5) may be solved for the correspond- 
Thus, we obtain 


nre(1 + cos 6) 
ry + re — c(sec? ¢ + tan? ¢) 


ing major axis 2a,,*. 


(2.2.7) 


2a," = 


Finally, the energy change AE,,* needed for this tra- 
jectory is found by substitution into Eq. (2.1.7). We 
have, after some simplification, 


AE,,* = 

c[(sec ¢ — tan — 2 sec ¢ tan 
+ cos 

(finite a,,*) 


(2.2.8) 


When AE+ is a monotonically decreasing tunction of 
a, the minimum value occurs for infinite a. 


The tra- 


jectory is parabolic and 
AE,,* = 3 — (2re/c) sin (Vr, (n + en c)+ 
(infinite a,*) (2.2.9) 


For the special case in which 6 = 7, the trajectory is 
of the Hohmann type and the minimum value of AF, is 
AExg, as obtained from Eq. (2.2.1). 

Fig. 2-3 gives plots of (Vre/Ve)»,* = V AE,,* as 
a function of 6for a voyage from Earth to each of the 
other planets of the solar system. The curves for the 
two most remote planets, Neptune and Pluto, are not 
shown because, with the scale used, they would be in- 
distinguishable from the curve for Uranus. The sec- 
tions of the curves for Jupiter, Saturn, and Uranus, 
corresponding to parabolic trajectories as the minimum- 
velocity paths, are characterized by broken lines. 


(2.3) Nonstop Round-Trip Interplanetary Trajectories 


Consider, as a specific mission for an interplanetary 
voyage, placing a spaceship in an orbit that passes 
within a few thousand miles of another planet and sub- 
sequently returns to Earth. The problem of deter- 
mining a suitable one-way trajectory may be solved 
with relative ease, using the material thus far devel- 
oped. The added complication of requiring the ve- 
hicle to return to Earth without additional propulsion 
(except that needed to correct for navigational inaccu- 
racies) would not contribute significantly to the diffi- 
culty of obtaining a solution were it not for the de- 
flection of the orbit caused by the gravitational field 
of the planet as the spaceship passes. Before treating 
the effect of perturbations introduced by the destination 
planet, let us consider first the nonstop round-trip 
trajectory, with the sun as the only gravitational force 
affecting the orbit. 

The simplest possible round-trip trajectory would 
be an orbit whose period is a multiple of the Earth’s 


period. For purposes of illustration and to obtain an 
URANUS 
SATURN 
ELLIPTICAL TRAJECTORY 
14+ — PARABOLIC TRAJECTORY 
10 
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Fic. 2-3. Minimum departure velocity as a function of polar 
angle 6. 
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Fic. 2-4. Earth to Mars departure velocity vs. semi-major axis 
of transfer ellipse for various values of the polar angle 6. 


appreciation for the magnitude of the velocities in- 
volved, let us investigate the possibilities of a spaceship 
orbit with a period of one year, that intersects the orbits 
of both Earth and Mars. From Eq. (1.6.1) we 
see that the semi-major axis a of such an orbit must be 
the same as that for the Earth—i.e., one astronomical 
unit (A.U.). Since the smallest value that a can 
assume is that for minimum energy, we find, from Eq. 
(1.2.2), that the largest possible linear distance be- 
tween the point of departure and the first crossing of 
the Martian orbit is 1.48 A.U., assuming the radius of 
the Martian orbit to be 1.52 A.U. Therefore, the 
spaceship must reach the Martian orbit after travers- 
ing an angular distance @ of not more than about 68.3°. 
From Eq. (2.1.7), we see that, for a one-year orbit and 
for the maximum possible value of 6, the required de- 
parture velocity is approximately 61,000 ft. per sec. 
Orbits with a one-year period requiring less velocity do 
exist for smaller values of 0; however, from the curves 
shown in Fig. 2-3, we see at once that this possibility 
ceases for angles less than about 35°. For example, we 
need a velocity of some 53,600 ft. per sec. for 6 = 58°; 
but, for 6 = 48°, the required velocity for a one-year 
orbit is more than 95,400 ft. per sec. 

The velocity requirements for an orbit whose period 
is 1.5 years are considerably more relaxed. In this 
case, the spaceship makes two circuits of the sun while 
the Earth makes three. The necessary semi-major 
axis is about 1.31 A.U. and calculations similar to 
the preceding show that such an orbit is possible for 
any value of 6. In Fig. 2-4, the lower portions of the 
Earth-to-Mars departure-velocity curves, as a function 
of the semi-major axis of the associated ellipse, are 
plotted for various values of 8. From these curves, we 
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see explicitly how the velocity requirements for a 
1.5-year orbit depend upon the polar angle 6. The 
minimum departure velocity for such a trajectory is 
approximately 11,200 ft. per sec. and the corresponding 
value for @ is about 142°. Fig. 2-5 gives the time of 
flight from Earth to Mars as a function of the semi- 
major axis, plotted for several values of 6. From these 
curves, the time to reach Mars on the minimum-depar- 
ture-velocity trajectory having a 1.5-year period is seen 
to be approximately 185 days. 

It is, of course, possible to have round-trip orbits 
that return to Earth after a non-integral number of 
years. However, we shall not pursue this point further 
before considering the additional effect of disturbances 
induced by the destination planet. 

The perturbations experienced by the orbit of a 
spaceship, when in the vicinity of a planet, depend upon 
the relative velocity with which the vehicle overtakes 
or is overtaken by the planet and the distance separat- 
ing the two at the point of closest approach. In the 
absence of any gravitational fields other than the 
planet’s own, the spaceship would approach the planet 
along a hyperbolic path. At a sufficiently great dis- 
tance the motion would be essentially along the asymp- 
tote and would have a velocity J. relative to the planet 
given by 


= 


according to Eq. (1.6.3). Here, we have denoted the 
semi-major axis of the hyperbola by a, and the gravi- 
tational constant of the planet by up. 

Referring to Fig. 2-6, we define 6 to be the angle be- 
tween the asymptote and the conjugate axis of the 
hyperbolic path of approach, «, to be the eccentricity 
of the orbit, D the distance between the vertex and the 
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Fic. 2-5. Earth to Mars time of flight vs. semi-major axis of 
transfer ellipse for various values of the polar angle @. 
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focus, and d the distance between the vertex and the 
suriace of the planet. The vertex is, of course, the 
point of closest approach of the spaceship to the planet, 
and we have the relationship 


D = — 1) = — 1) 
Solving for ¢, and noting that e«, = csc 6 we obtain 
sin 6 = + (DV."/up)| (2.3.1) 


Therefore, the total effect on the velocity of the space- 
ship, after contact with a planet, is simply a rotation 
in the plane of motion of the relative component of 
velocity by an amount 26. The direction of rotation 
of the relative velocity vector will tend to increase or 
decrease the absolute velocity according to whether 
the spaceship passes behind or ahead of the planet, 
respectively. 

The mass ratio of the sun to Mars is 3,093,500. 
Therefore, the attraction of the planet does not become 
as large as one hundredth of the attraction of the Sun 
until the spaceship is within about 800,000 miles of 
Mars. By a direct calculation, using Eqs. (1.5.9), 
(2.1.2), and (2.1.3), we find that a vehicle moving in a 
1.5-year orbit that intersects the Martian orbit at 
@ = 140° has a velocity component in the direction of 
the motion of Mars of about 48,800 m.p.h. and a radial 
component away from the sun of approximately 8,770 
m.p.h. Since the orbital velocity of Mars is about 
54,000 m.p.h., the spaceship is being overtaken by the 
planet and the relative velocity between the two 
bodies is more than 10,000 miles per hour. Thus, the 
vehicle is within the critical 800,000-mile distance ot 
the planet for less than 6.5 days. For other round-trip 
orbits, the relative velocity is higher, so that the 
period of time during which Mars can influence the orbit 
will be correspondingly smaller. 

Because of the above considerations, it seems rea- 
sonable, as a good approximation to the true state of 
affairs, to view the effect of Mars on the trajectory as 
simply an impulse in velocity applied at the instant 
the spaceship crosses the Martian orbit. In computing 
the magnitude and direction of the velocity impulse, we 
shall use Eq. (2.3.1) to obtain the turn angle 26, 
with |”,, taken as the relative velocity of the spaceship 
with respect to Mars, determined at the point of inter- 
section of the two orbits. 

Considering again the problem of a vehicle moving 
in a 1.5-year orbit that intersects the orbit of Mars at 
6 = 140°, let us assume that the closest approach to 
the surface of the planet is 3,000 miles and that the 
ship passes ahead of Mars. Then, following the pro- 
cedure outlined, the relative velocity is 10,245 m.p.h. 
and the turn angle 26 is 23°. Therefore, the effect 
of the contact is to reduce the absolute velocity from 
49,540 to 46,150 m.p.h. Furthermore, the period of 
the orbit is reduced from 1.5 to 1.31 years, as deter- 
mined from Eqs. (1.6.2) and (1.6.1). These are, of 
course, sizeable changes, which must be taken into 
account in calculating round-trip interplanetary re- 
connaissance trajectories. 
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Fic. 2-6. Motion of spaceship in vicinity of Mars. 


(2.4) Reconnaissance Trajectories for the Planet Mars 


Using our simple model of the solar system, we shall 
now describe a method for computing nonstop tra- 
jectories from Earth to a planet and return. For 
definiteness, we shall choose the planet Mars as the 
target and consider, as allowable spaceship orbits, 
elliptical trajectories requiring reasonable departure 
velocities. The calculation technique to be described 
is not necessarily the most efficient considering the 
simplicity of our model; however, the ideas are readily 
adapted to the more difficult three-dimensional problem 
that will be treated in Section (3). 

Let Vz be the orbital velocity of the Earth and let a 
departure velocity magnitude Vr, from the Earth's 
orbit be specified. We now choose a value 6p for the 
heliocentric angle @ between the point of departure on 
the Earth’s orbit and the point on the Martian orbit at 
which we intend to intercept Mars. Then, with AE = 
(Vre/ Ve)? and 6 = 6p, Eq. (2.1.7)f is solved for the 
corresponding value of the semi-major axis dp of the 
departure orbit. Obviously, a solution is not always 
possible. On the other hand, the equation may have 
two roots, each of which corresponds to a satisfactory 
solution to the problem of determining a departure 
orbit. 

The next step in the procedure is to compute the 
time 7 required for the journey from Earth to Mars. 
For this purpose, we may use the following convenient 
formula, which is universally valid if @ is not an integral 
multiple of x. Thus, if we define N to be the integral 

7 Actually, Eq. (2.1.7) holds only for positive values of sin @. 
Therefore, if 6 lies in the third or fourth quadrant, we must use 
the symmetry inherent in our simple model—i.e., for any in- 
teger k, we observe that @ and 2kx — 6 require identical departure 
orbits for the same velocity. 
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part of 0/27, we may use Eqs. (1.6.6) and (1.6.7), described, we will either find a return trajectory that 
together with the properties of symmetry, to show that contacts the Earth or conclude that none exists. 
, Assuming that a return trajectory has been found, 
T = (P/2) {2N +1 + [sgn (sin 0)/x] X 
[+( a pe is — de 8)]} (2.4.1) it remains to determine whether the gravitational 
” field of Mars is sufficient to accomplish the transfer 
The positive or negative sign is used, respectively, of the spaceship from the departure orbit to the return 
according to whether the root a@p was obtained from orbit. Using Eqs. (2.1.2) and (2.1.3), we may compute 
the AE, branch of the curve or the AE_ branch. The the necessary angle 26 through which the relative 
function sgn is defined as +1, 0, or —1, according to velocity vector must be rotated in order to change 
whether the argument is, respectively, positive, zero, orbits. Then if Vry is the magnitude of the relative 
or negative. velocity vector and wy is the gravitational constant 
With the time for the trip from Earth to Mars deter- of Mars, we may compute the required passing distance 
mined, the location of Mars at the time of departure D from 
is fixed. Thus, the orientations of the planets, and ee iy 
hence the time at which the voyage may start, is a con- D = (um/Vem’*)(ese 6 — 1) (2.4.3) e 
sequence of the selected departure orbit. ae If D has a reasonable value, considering the objectives 
The relative velocity Vey of the spaceship with of the reconnaissance mission, we will have a solution (2 
respect to Mars can be calculated at the point of inter- to the round-trip trajectory problem. : 
section with the Martian orbit from the relationship A digital computer program has been prepared that 
Van? = Vu — — 2V (2.4.2) mechanizes the preceding computational procedure. (6 
By specifying the departure velocity and an angular 
where Vy is the orbital velocity of Mars and /p is the increment used to generate values of 6p, the computer ce 
semi-latus rectum of the departure orbit. will systematically find all of the corresponding round- re 
Setting aside for the moment the effect of Mars on trip trajectories. For purposes of illustration, the re- s) 
the orbit, let us consider next the problem of returning sults of this program are summarized in Table 1 for a s) 
to Earth. The procedure for determining an appro- departure velocity equal to 0.13 times the Earth’s pt 
priate return trajectory that intersects the orbit of the orbital speed and an angular increment of 10°. Cer- n 
Earth is the same as used in obtaining the departure tain arbitrary bounds were placed on the program to fi 
trajectory. The magnitude of the relative velocity with prevent the computer from generating a trajectory re- 
which the ship leaves Mars is the same as that with quiring more than 3.5 years for the round trip or passing fie 
which it arrived, regardless of the perturbations in- closer than 500 miles to, or farther than 30,000 miles pe 
duced by the contact. Therefore, we may proceed as from, the surface of Mars. The notation and units ck 
before by selecting the heliocentric angle 6g measured used in Table 1 are listed as follows: ay 
from the intercept point on the Martian orbit to the ti 
point on the Earth’s orbit where we may hope to re- 6) = heliocentric angle through which the space- of 
establish contact with the Earth. Again, we use Eq. ship moves from Earth to Mars, deg. a 
(2.1.7), with AE = (Vru/ Vy)? and with 7; and rz now ap = semi-major axis of the departure ellipse, th 
taken as the radii of the orbits of Mars and Earth, A. U. i, 
respectively, to solve for the semi-major axis dz of the Ep = eccentricity of the departure ellipse 
return orbit. If a solution exists, Eq. (2.4.1) may be Tem = time for trip from Earth to Mars, days or 
used as before to determine the time for the return trip. Vew = relative velocity of vehicle with respect to nc 
With the time for the complete round trip known, Mars at point of intersection with the M 
it is a simple matter to determine the location of the Martian orbit, ft. per sec. ve 
Earth at the instant the spaceship returns to the orbit 26 = angle of rotation of the relauve velocity ce 
of the Earth. By systematically varying the choice vector resulting from contact with Mars, be 
of the angle 6r and repeating the calculations just deg. pl 
Trip 8p, ap, Tem, Vem, 26, d, Or, ar, Tue, diem, Tror, ay 
No. deg. A.U. €D days _ft./sec. deg. miles deg. A.U. €R days deg years pc 
1 130 1.379 0.275 164 20,300 85 6,600 528 1.297 0.301 869 44 2.828 M 
2 136 6.7 9,100 632 1.314 0.296 974 3.115 sh 
3 230 “ S 428 sia 17.3 1,800 434 1.228 0.333 612 6 2.846 
4 230 11.3 4,300 563 1.274 0.311 742 3.203 A: 
5 490 ~ a 756 i 14.3 2,800 283 1.250 0.322 393 — 266 3.147 . 
6 140 1.359 0.266 175 19,100 6.5 11,100 534 1.301 0.286 874 48 2.872 Ue, 
7 140 5.0 15,300 625 1.314 0.281 966 3.124 di 
8 220 i fa 403 % 13.1 4,000 442 1.251 0.308 634 9 2.838 
9 220 “¢ 8.5 7,800 555 1.285 0.293 748 3.152 
10 500 > ss 754 is 19.9 1,700 130 1.207 0.332 250 —255 2.750 ax 
11 500 10.4 5,800 275 1.271 0.299 397 3.152 E 
12 210 1.334 0.255 375 17,400 7.5 11,500 455 1.278 0.278 665 14 2.846 
13 510 $ ig 751 iy 9.8 8,100 157 1.262 0.285 291 —243 2.853 fo 
14 510 5.0 18,600 260 1.295 0.270 396 3.140 Q 
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d = minimum passing distance from the surface 
of Mars, miles 

Or = heliocentric angle through which the space- 
ship moves from Mars to Earth, deg. 

de = semi-major axis of return ellipse, A.U. 

ER = eccentricity of return ellipse 

Tue = time of trip from Mars to Earth, days 

Agm = heliocentric angle between Earth and Mars 


at instant of departure, in degrees meas- 
ured in the counterclockwise direction 
from Earth to Mars 

Tror = total time for the round trip, years 


In Section (3), certain of these trajectories will be 
compared with their three-dimensional counterparts. 


(3) Planetary Reconnaissance Trajectories in a 
Three-Dimensional Model of the Solar System 


(3.1) Coordinate Systems, Angles and Directions 


As a preliminary step in the development of a pro- 
cedure for determining nonstop, round-trip planetary 
reconnaissance trajectories in a true model of the solar 
system, it will be convenient to introduce the various 
systems of coordinates that will be needed. Wherever 
possible we shall adopt the coordinate systems and 
notation that have been universally accepted in the 
field of celestial mechanics. 

The position of the Earth at any time will be speci- 
fied in terms of a set of heliocentric rectangular com- 
ponents. As shown in Fig. 3-1, the x and y axes are 
chosen in the plane of the ecliptic, with the positive x 
axis in the direction of the vernal equinox. The posi- 
tive y axis is in the general direction of the perihelion 
of the Earth’s orbit, and the z axis is chosen to make 
a right-handed coordinate system. Unit vectors in 
these three directions will be defined, respectively, as 
i,, i,, and i,. 

The line of intersection of the plane of the Martian 
orbit with the plane of the ecliptic is called the line of 
nodes. The ascending node AN is the point at which 
Mars crosses the ecliptic with a positive component of 
velocity in the z direction. The longitude of the as- 
cending node as measured from the vernal equinox will 
be denoted by 2. The angle of inclination of the orbital 
plane of Mars to the ecliptic is called 7. 

To specify the location of Mars, a different set of 
heliocentric coordinates will be used. The £& and n 
axes are selected in the Martian orbital plane with the 
positive ¢ axis in the direction of the perihelion of the 
Martian orbit. The y and ¢ axes are then chosen, as 
shown in Fig. 3-1, to make the system right-handed. 
Associated with the three axes are the unit vectors 
iz, ip, and i;. The & axis makes an angle w with the 
direction of the ascending node. 

Let dz, €gz, and /g be, respectively, the semi-major 
axis, the eccentricity, and the semi-latus rectum of the 
Earth’s orbit, and correspondingly define ay, €4, li 
for Mars. Furthermore, denote by IIy and ly = 
Q + w the respective longitudes of the perihelions of 
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Fic. 3-1. Fixed heliocentric rectangular coordinates for Earth 
and Mars orbits. 


the Earth and Martian orbits measured with respect 
to the vernal equinox. Finally, if =, and Zy are the 
longitudes of the planets at the epoch, then, correspond- 
ingly, at a time 7 years from epoch their mean anoma- 
lies, My and My are given by 


Mer = 2eT + =E 
Mu (24T Px) + Zu — 


(3.1.1) 
(3.1.2) 
where Py, the period of Mars in years, is given by 
Py = (3.1.3) 
From the mean anomalies, we may compute the 
eccentric anomalies Ey and Ey, by means of the Fourier- 
Bessel expansions 


Ex = Mg (1/n)J,(neg) sin (3.1.4) 


n=1 
Ey = My +2> (1/n)J,(neu) sin nMy (3.1.5) 
n=1 
where J, is a Bessel function of the first kind of order n. 


The position vectors rg and ry at time 7 may then be 
determined from 


fe = [ag cos IIg(cos Ex — ex) — 
V arly sin II, sin Eg li, + 
[az sin IIz (cos Eg — €g) + 
V cos IIg sin Egli, (3.1.6) 


= ay(cos Ey — + Vauly sin Eyi, (3.1.7) 


In order to obtain the components of the Martian 
position vector ry in the ecliptic coordinate system, one 
may simply pre-multiply the vector ry, regarded as a 
column vector, by the matrix 


h bh ds 


=| mM. Ms 
Nn Nz 


(3.1.8) 


cos 2 cos w — sin 2 sin w cos 7 


where /; 


l, = —cosQ2 sin w — sin 2 cos w cost 
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< OF MARS 


/ ORBIT 
OF EARTH 


\ ~ ANGLE OF INCLINATION OF TRAJECTORY PLANE TO PLANE OF THE ECLIPTIC 
ANGLE OF INCLINATION OF TRAJECTORY PLANE TO PLANE OF THE MARTIAN ORBIT 


Fic. 3-2. Local polar coordinate systems of Earth, Mars, and 
spaceship orbits. 
= sinQ sini 
m, = sin Q cos w + cos 2 sin w cos 7 
m, = —sin Q sin w + cos 2 cos w cos 7 
m3; = —cosQ sini 
n = sinw sin7z 
= cosw 
nN; = cost 


Consider now the problem of connecting a point P 
on the Earth's orbit and a point Q on the Martian orbit 
by an ellipse whose focus F is sun-centered, as illus- 
trated in Fig. 3-2. Let fy, and fy, represent the vector 
positions of Pand Q. Then the cosine of the angle 6 
between fr, and fy, is obtained from 


cos = (3.1.9) 


Since @ is measured in the direction of motion of the 
planets, the sine of the angle, with appropriate con- 
sideration for sign, is given by 


sin = sgn (rz X V1 — cos? 6 (3.1.10) 


The unit vector i,, normal to the plane of the ellipse and 
having a positive component in the é, direction may then 
ke obtained from 


i, = (re x fu) sin 0 (3-5. 1 1) 


When computing the vector velocities of the Earth 
and spaceship at the point P, it will be convenient to 
use a local polar coordinate system centered at P, 
with unit vector i, in the positive direction of ry, and 
with i, chosen in the plane of the ecliptic and such that 
i,, i, and i, form the right-handed triad shown in 
Fig. 3-2. The plane of the ellipse is inclined at an 
angle x to the ecliptic, with the angle positive in the 
direction of a right-handed rotation about i,. Hence, 


cos xX = ini, (3.1.12) 


One may readily see that x is positive if, simultane- 
ously, @ is less than 180° and the point Q is above the 
ecliptic or, alternately, @ is greater than 180° and Q is 


These 
facts may be summarized conveniently by means of 
the equation 


below the ecliptic; otherwise, x is negative. 


sin x = sgn (f,,-i, sin 8) Vi — cos? x (3.1.13) 


Similarly, at the point Q we introduce a set of unit 
vectors i, and i,, as shown in Fig. 3-2. The angle » 
is the inclination angle of the trajectory plane with 
respect to the plane of the Martian orbit and is positive 
in the direction of a right-handed rotation about i,. 
Therefore, 


cos v = i, (3.1.14) 


Again one sees that v is negative if, simultaneously, 6 is 
less than 180° and the point P is above the Martian 
plane or, alternately, @ is greater than 180° and P is 
below the Martian plane; otherwise, vy is positive. 
Hence, 


sin y = sgn sin 6) V1 — cos? (3.1.15) 


We shall close this preliminary section with the der- 
ivation of a criterion for determining the sign of the 
radial component of the vector velocity of the spaceship 
at the points P and Q. For this purpose, let us define 
the vector i; as a unit vector from the center of force 
F and in the direction of the vacant focus of the ellipse. 
Since i, lies in the plane containing both ry and ry, it 
may be expressed as a linear combination of these 


two vectors—1.e., 
+ Fyfy (3.1.16) 


By taking the scalar product of i, with ry and fy, in 
turn, we have 

= Ferg’ cos 6 

= Ferery cos 0 + Fyry? 


But if % is the angle between fr, and i,, and / and ¢ are, 
respectively, the semi-latus rectum and eccentricity 
of the ellipse, then, using the polar-coordinate equa- 
tion for an ellipse, given by Eq. (1.5.1), 


ip-f, = COS = (rz — 
iv fu = COS = ), € 
Hence, solving for Fy and Fy,, we obtain 


Fy = (1/erg sin? — (//rg)] — 


[1 — (//ryz)] cos (3.1.17) 
Fu => (1/eray sin? { [1 | 
[1 — cos 6{ (3.1.18) 


At the point P, the radial component of the vector 
velocity of the spaceship and the component of the 
vector fz X fr, in the i, direction have identical signs. 
Therefore, it follows from Eqs. (3.1.16) and (3.1.11) 
that the spaceship will be moving away from or toward 
the sun according to whether the sign of Fy sin @ is 
respectively positive or negative. Similarly, at the 


point Q, the corresponding criterion is the sign of 
—F, sin 6. 
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(3.2) Procedure for Calculating Round-Trip Trajectories 


Essentially the same technique described in Section 
(2.4) for determining round-trip reconnaissance tra- 
jectories in the simplified model of the solar system 
may be used in the three-dimensional version. One 
fundamental difference occurs in the method of esti- 
mating the point on the Martian orbit where planetary 
contact can occur and, subsequently, the point at 
which the spaceship will return to the Earth. With 
circular symmetry no longer a characteristic of the 
model, we are not free to choose arbitrarily the helio- 
centric angle 6. However, if we postulate a departure 
time, 7,, measured in years from the epoch, we may, 
instead, estimate the time 7’, required for the voyage 
from Earth to Mars and the time 7, for return. 
The net effect is, of course, the same, in that the angles 
6p and 6g are direct functions of these times. As an 
aid in the estimation of the times Tp, Ty, and T yz, 
the results obtained for the two-dimensional model can 
be used to advantage. 

The two methods also differ in the manner of deter- 
mining the semi-major axis of the spaceship orbit. We 
bave no three-dimensional analog of Eq. (2.1.7); in- 
stead, since the times of flight are fixed by our initial 
estimate, the time-of-flight equation (2.4.1) must serve 
as the means for calculating ap and ap. 

If, in addition to the departure time, 7p, we specify 
the relative departure velocity magnitude, Vz, it will 
be necessary to repeat the calculation of ap while 
systematically varying the time 7 in order to ob- 
tain a trajectory that satisfies this additional require- 
ment. Under any circumstances, this procedure will 
be mandatory for calculating the return-trip trajectory 
since the vehicle must leave Mars with the same rela- 
tive velocity magnitude, V.,, with which it arrived. 

Before detailing the step-by-step process, let us turn 
our attention briefly to the problem of solving the 
time-of-flight formula for the semi-major axis a. 
Curves of the time of flight, 7, as a function of a were 
computed trom Eq. (2.4.1) and are plotted in Fig. 3-3 
for several values of 6. If the angle 6 is less than 360° 
(N = 0), Eq. (2.4.1) defines a as a single-valued func- 
tion of 7. Therefore, if a trajectory connecting the 
points P and Q is possible for a given value of 7, it is 
unique. On the other hand, if 6 is greater than 360°, a 
is a doubled-valued function of 7. Thus, correspond- 
ing to each value of 7 that is sufficiently large to ensure 
a solution, two different trajectories are obtained. 

For the case in which NV = 0, the procedure for ob- 
taining a is straightforward. We first compute the 
time of flight 7,, for the minimum-energy path, for 
which a = ad» or, equivalently, a = 7m, from the formula 


Tn = (Pm/2){1 — [sgn (sin 0)/2](Bm — sin 
(3.2.1) 


Then by comparing the given value of 7 with T,,,, we 
find that the upper or lower sign is appropriate in Eq. 
2.4.1) according to whether 7 is, respectively, less 
than or greater than 7,,. However, if 7 is less than 


T (yr) 


(AU) 


Fic. 3-3. Time of flight from Earth to Mars vs. semi-major axis 
of transfer ellipse for various values of the polar angle @. 


the time 7, required to travel the parabolic path 
from P to Q, where 


T, = (1/3) V2 ae [s*/? — sgn (sin @)(s — c)*/?] 


then no solution for elliptical trajectories is possible. 
With the feasibility of the elliptical orbit confirmed and 
the appropriate sign determined for Eq. (2.4.1), a 
simple iteration will produce the unique root. 

The situation is somewhat more complex when JN is 
not zero. If a relevant two-dimensional solution to 
the trajectory problem is available, it can be used to 
advantage in selecting the proper branch of the curve 
a asa function of 7. In the absence of any such infor- 
mation, one can always proceed with the computation, 
using both roots, and reject whichever proves unsatis- 
factory. 

Returning now to the main problem, we shall out- 
line a procedure for calculating a nonstop trajectory 
for a voyage from Earth to Mars and return, to begin 
on a specific date with a prescribed relative departure 
velocity magnitude. 

(1) A departure time, 7p, is specified that is meas- 
ured in years from the epoch. 

(2) The mean anomaly of the Earth at the point of 
departure, M,(Tp), is computed from Eq. (3.1.1). 
Then, from Eg. (3.1.4), the corresponding eccentric 
anomaly of the Earth, Eg(Tp), is obtained. Finally, 
the vector position of the Earth at the time of depar- 
ture, f¢(Tp), is computed from Eq. (3.1.6). 

(3) The vector velocity of the Earth at time Tp is 
obtained from 


V;,(Tp) = sgn [sin E,(Tp)| x 


V — — (1/ae)} + 
(3.2.3) 
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(4) An estimate of the time 7, in years required 
for the voyage from Earth to Mars is made and the 
mean anomaly of Mars, My(Tp> + Tx), for the point 
at which contact with the spaceship will occur is com- 
puted from Eq. (3.1.2). The corresponding eccentric 
anomaly, Ey(Tp + Ty), and vector position, p 
+ Tym), are then obtained from Eqs. (3.1.5) and (3.1.7). 
(5) The heliocentric angle @p, through which the 
spaceship moves on the departure orbit, is calculated 
using Eqs. (3.1.9) and (3.1.10), from which the linear 
distance c, measured from the point of departure P to 
the point of arrival Q at Mars, is obtained from 


c? = re? + ry? — 2rery cos Op (3.2.4) 


(6) From the time-of-flight equation (2.4.1), with T = 
T zm, the semi-major axis @p of the departure orbit is 
obtained. Then the semi-latus rectum /p is computed 
from Eq. (1.5.9), with the choice of upper or lower sign 
made to correspond with the sign used in Eq. (2.4.1) 
to produce the root ap. 

(7) The vector i,, normal to the plane of the departure 
trajectory is computed from Eq. (3.1.11) and the result 
used in Eqs. (3.1.12) and (3.1.13) to obtain the inclin- 
ation angle x. 

(8) The velocity vector Vpy of the spaceship at the 
point of departure is computed from 


Voge = sgn (Fy sin 0p) X 
V [(2rz — lp)/rz*| — (1/ap)}i, + 
(WV (cos xip + sin xi) (3.2.5) 
where Fy is defined in Eq. (3.1.18). 
(9) The vector velocity Vez of the spaceship rela- 


tive to the Earth at the point of departure P may then 
be computed as the vector difference 


Vee = — 


(10) The magnitude of the relative velocity Vez may 
now be compared with the prescribed value. If the 
two differ by more than a tolerable amount, a new esti- 
mate for 7,4 must be made and steps (4) through (9) 
repeated. 

(11) When an Earth-to-Mars orbit has been found 
that satisfies the initial departure conditions, we may 
then determine the relative velocity of the spaceship at 
Q, the point of contact with Mars. For this purpose, 
the vector velocity of Mars at time Tp + Ty, is com- 
puted from 


VulTpd + Tem) = sgn [sin Ey(Tp + Tem)] X 
V us{ — — (1/au)}i, + 
(3.2.7) 


and the vector velocity of the spaceship at the point 
Q from 


Vom = sgn sin 6p) X 
V ust — — (1/ap)}i, + 
(VW (cos vi, + sin vi,) (3.2.8) 


where Fy is defined in Eq. (3.1.17) and the inclination 


(3.2.6) 
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angle v is obtained using Eqs. (3.1.14) and (3.1.15). 
Finally, the relative velocity Vey is computed as the 
vector difference 


Vew = Vow — Vul(To + Tem) (3.2.9) 


(12) To obtain a return trajectorv, the procedure is 
precisely the same. The spaceship departs from Mars 
at time Tp + Tyg and arrives at Earth Ty, years 
later. As before, the time for the trip from Mars to 
Earth is a result of the iteration process used for de- 
termining a trajectory whose relative velocity magni- 
tude Vey at Mars is prescribed. 

(13) As a final step, we must determine the orien- 
tation of the spaceship with respect to Mars that must 
be effected during the period of contact in order to 
accomplish the transfer from the departure orbit to the 
return orbit. Let Vryp and Vrwr be, respectively, 
the relative velocities with respect to Mars of the ve- 
hicle on the departure orbit and the return orbit. The 
turn angle 26 may be computed from 


| x | 


Vem’ 


sin 26 = (3.2.10) 
and the passing distance D from Eq. (2.4.3). 

If the vector product in Eq. (3.2.10) has a positive 
component in the ¢ direction, we shall say that the 
spaceship passes ahead of Mars; otherwise, the ship 
passes behind the planet. Furthermore, this vector 
product defines the orientation of the plane of relative 
approach to Mars. 


(3.3) Comparison of Results With the Simplified Model 


The computational procedure for determining three- 
dimensional round-trip trajectories has been programed 
for a digital computer. The input to the program con- 
sists of a departure time 7p and a departure velocity 
Vere, together with estimates of the times required for 
the trips to and from Mars. 

Several of the two-dimensional trajectories found in 
Section (2) were recalculated in the three-dimensional 
model. The same departure velocity of 0.13 times the 
Earth’s orbital speed, or 12,709 ft. per sec., was used. 
The time of departure Tp was selected so that the dif- 
ference of the mean longitudes of the planets would be 
the same as the heliocentric angular difference A gy be- 
tween Earth and Mars tound in the two-dimensional 
solution. The time-of-flight estimates were likewise 
taken directly from the results in the simplified model. 

For comparison purposes, the corresponding tra- 
jectories are represented schematically in Figs. 3-4 
through 3-9. In each case, the circles shown are sun- 
centered, with radii equal to the mean distances of the 
Earth and Mars. The major axes of the true elliptical 

orbits of the planets are shown as broken lines, with the 
perihelions and aphelions appropriately marked Pz, 
Ary, and Py, Ay. The line of nodes is also shown as 


a broken line, with the ascending and descending nodes 
labeled AN and DN. The positions of the planets at 
the times of departure, arrival at Mars, and return to 
the Earth are denoted by the letters E and M with 
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respective subscripts D, A, and R. The point of de- 
parture from the Earth E> is approximately the same 
for two models; however, the points of arrival at Mars 
and return to Earth are different. The primed symbols 
correspond to the three-dimensional model, while the 
unprimed symbols correspond to the two-dimensional 
model. The major axes of the departure and return 
orbits for both cases are shown as solid lines, with the 
perihelion and aphelion points labeled and the positions 
of the vacant foci marked. 

Together with the diagrams are shown some of the 
data for the three-dimensional trajectories. These 
data may be compared with Table 1 of Section (2). 
As a result of this comparison, certain general con- 
clusions may be drawn, which are discussed below. 


The most striking difference between the two models 
is seen for those trajectories in which the point of de- 
parture from the Earth and the point of arrival at the 
planet are nearly 180° apart. For the case of co-planar 
planetary orbits we have seen that this configuration 
of the planets gives rise to the absolute minimum- 
energy transfer ellipse. In the three-dimensional 
model, even though the angle between the two planes 
is less than 2°, the situation is considerably altered. 
If the diametrical position of the points of departure 
and arrival lie on the line of nodes, the orbital plane of 
the spaceship may be chosen arbitrarily. On the other 
hand, since this plane must always contain the three 
points determined by the positions of Earth, sun, and 
Mars, a large angle of inclination to the ecliptic can re- 
sult when these three points are not distributed on the 
line of nodes. In particular, when the diametrical 
position of the planets is at right angles to the line of 
intersection of the orbital planes of the planets, the 
trajectory plane will be at right angles to the plane of 
the ecliptic. 

This point is rather dramatically illustrated when 
trajectories No. 1 and No. 6 are contrasted. In many 
respects the two seem alike. For each, the spaceship 
leaves at roughly the same time, contacts Mars at 
essentially the same point in space after some six 
months, and each returns to Earth at approximately the 
same position. However, we note that for the return 
orbit the points of arrival at Mars and return to Earth 
are 184.76° apart for trajectory No. 1 and 182.75° 
apart for trajectory No. 6. The inclination angles of 
the return orbits are, respectively, 4.671° and 14.481°. 
We further note, in the latter case, that a 44° rotation 
of the relative velocity vector at Mars must be effected 
in order to transfer from the departure orbit to the 
return orbit. The gravitational field of Mars is inade- 
quate to accomplish this transfer, as seen from the fact 
that a mathematical solution is obtained only by per- 
mitting the spaceship to pass more than 1,200 miles 
beneath the Martian surface. 

Still another difference between the two models is 
apparent in the manner of approaching Mars. In the 
two-dimensional model, all motion takes place in a 
single plane. Even in the three-dimensional model, 
except for the singular situation noted, the inclination 


DEPARTURE TIME: DECEMBER 13, 1964 
RETURN ORBIT: 


DEPARTURE ORBIT: 


%® = 
Veep - 
Vemp 
VreD * 
Vaup 
te (Tp) = 

ty(Tp+ Tey) = 

Tp = 

Tew 

€p - 

xp * 


135.79° 

69437, + 10591ig+ 1065i, 
17016i, - 7597 + 1684i- 
12709 ft/sec 

18711 ft/sec 

0.9843 A.U. 

1.5800 A.U. 

6.951 years 

0.5262 years 


1.3091 A.U. 
0.2543 


0.555° 
1,409° 


Fic. 3-4. 


= 544.76° 
Vaer ~ ~7063i,+ 10128i 4+ 8882, 


Vawr 167211, 767914, - 33917 
Vrer = 15210 ft/sec 
Tye = 2.3669 yeors 
Op = 1.3076 AU. 
= 0.2532 
xp = 4.6719 
vp = 2.840° 
PLANETARY CONTACT: 
25 = 15.619° 
d = 3110 miles 


THE VEHICLE PASSES AHEAD OF MARS IN A 
PLANE INCLINED 87.266° TO THE MARTIAN 


ORBIT. 


Trajectory No. 1. 


DEPARTURE TIME: 


DEPARTURE ORBIT: 


5p = 
Vrep 
RMD 
Vamp 
tg(Tp) = 
tylTp + Tew) 
Tp 
Tem 
Op = 
Xp * 


134.87° 

18653, + 124281 9+ 18897, 
226671, - 63387, + 12367, 
12709 ft/sec 
23569 ft/sec 

0.9854 A.U. 

1.6023 A.U. 

6.927 yeors 

0.4943 yeors 

1.3785 A.U. 

0.2853 

0.970° 

1.031¢ 


DECEMBER 4, 1964 
RETURN ORBIT: 


6p = 5$42.75° 
= - 74621, +7293ig+ 272667, 
Vaun 158737, ~ 8971, 149347 
* 29194 ft/sec 
Tyg = 2.3902 yeors 
og = 1.3185 AU. 
~ 0.2583 
xp = 14.481° 
- 12.7679 
PLANETARY CONTACT: 
25 = 44.204° 
d = —1252 miles 


THE VEHICLE PASSES AHEAD OF MARS IN A 
PLANE INCLINED 73.862° TO THE MARTIAN 


ORBIT. 
Fic. 3-5. Trajectory No. 6. 
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DEPARTURE ORBIT: 


6p = 134.87° 
18657, + 124287 4+ 18897, 
Vaup 226677, - 63387, + 12367, - 
Veep 12709 ft/sec 
Vamp = 23569 ft/sec 
te(Tp) = 0.9854 A.U. 
ty(Tp+ Tey) = 1.6023 AU. 
Tp = 6.927 yeors 
Trey = 0.4943 yeors 
Op = 1.3785 A.U. 
= 0.2853 
xp = 0.970° 
vp = 1.031¢ 
Fic. 3-6. 


Vrep 
Vamp 
VreD 
Vamp 

te (Tp) 
tulTp + Tew) 
Tp 
Tem 
€p 
XD 


DEPARTURE TIME: 


DECEMBER 4, 1964 
RETURN ORBIT: 
Op = 624.95° 
= 219471, + 8577ig+ 13017, 
21805;, - 87397, + 19167, 
23599 ft/sec 
Typ 26128 years 


Op = 1.3123 

= 0.3176 

= 0.690° 

vp = 1.656° 
PLANETARY CONTACT: 
25 = 6.420° 
d - 6611 miles 


THE VEHICLE PASSES AHEAD OF MARS IN A 
PLANE INCLINED 16.598° TO THE MARTIAN 
ORBIT. 


Trajectory No. 7. 


DEPARTURE TIME: FEBRUARY 28, 1965 


DEPARTURE ORBIT: 


219.91° 

61857, + 107497 4+ 27807, 
-24593i,, - 83687, 36281 
12709 ft/sec 

26230 ft/sec 

0.9906 A.U. 

1.4148 A.U. 

7.160 yeors 

1.0744 yeors 

1.3168 A.U. 

0.2565 

1.455° 

-2.713° 


Fic. 3-7. 


RETURN ORBIT: 
9p = 510.00° 
= 90541, + 102321 4+ 35431, 
Vawr , - 87303, - 4280; - 
Veer 14114 ft/sec 
Tye 1.9527 years 
ag = 1.3062 AU. 
0.2574 
xp = 1.868° 
vp = -3.215° 
PLANETARY CONTACT: 
25 = 1.706° 
d = 25519 miles 


THE VEHICLE PASSES BEHIND MARS IN A 
PLANE INCLINED -33.518° TO THE MARTIAN 
ORBIT. 


Trajectory No. 9. 


DEPARTURE TIME: APRIL 27, 1960 


DEPARTURE ORPIT: 


= 470.43° 
Vrep ~ 1796i,+ 12016 i 37291, 
Vaup 220477, 74097, + 613i, 
VrEp |2709 ft/sec 
Vamp = 23266 ft/sec 
te(Tp) = 1.0066 A.U. 


ty(Tp+ Tey) = 1.3828 A.U. 


Tp = 2.319 years 
Try ~ 1.9259 years 
= 1.3523 A.U. 


= 0.2556 
Xp = 1.957° 
vp = 0.442° 
Fic. 3-8. 


RETURN ORBIT: 

Oy = 287.34° 
172957, + 85777,- 35197, 
~ 210407, 98717, - 10927. 
Vrer ~ 19623 ft sec 
Typ = 1.1814 years 


Gp ~ 1.2755 A.U 

0.2634 

= 

vp = 0.812° 
PLANETARY CONTACT 
25 = 7.785° 
d = 5179 miles 


THE VEHICLE PASSES AHEAD OF MARS IN A 
PLANE INCLINED 32.082° TO THE MARTIAN 
ORBIT. 


Trajectory No. 11. 


DEPARTURE TIME: APRIL 1, 1960 


DEPARTURE ORBIT: 


Op = 478.71° 
~ - 35571, + 11535ig- 39767, 
Vamp 233751 77741 
Vrep = 12709 ft'sec 
Vamp ~ 24698 ft/sec 
te(Tp) = 0.9995 A.U. 


ry(Tp + Tey) = 1.3929 AU. 


Tp = 2.250 years 
Tey = 1.9235 yeors 


op, 1.3364 AU. 
€p 0.2527 
xp = -2.083° 
vp 1.293° 


RETURN ORBIT: 
= 201.61° 
~ -13801i, 93531, 
~ 216901,- 98171, + 65671 - 
VreR = 18859 ft/sec 
Typ = 9.9675 years 
Op = 1.2755 AU 
= 0.2512 
xr = 4.965° 
4.909¢ 
PLANETARY CONTACT 
25 = 12.736° 
d = 1666 miles 
THE VEHICLE PASSES AHEAD OF MARS IN A 
PLANE INCLINED 63.873° TO THE MARTIAN 
ORBIT. 


Fic. 3-9. Trajectory No. 13. 
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angles of the planes of motion of the spaceship are small. 
However. as viewed from the surface of Mars, the 
angular orientation of the plane of approach is cer- 
tainly not small. Thus in the vicinity of Mars, the 
true relative motion of the spaceship and planet is not 
at all adequately approximated by the two-dimensional 
model. 

A third point of contrast has to do not so much with 
the three-dimensional nature of the solar system as 
with the ellipticity of the planetary orbits. In the 
elementary model, the planets were assumed to move 
in concentric circles about the sun whose radii were 
taken as the corresponding mean distances. The ap- 
proximation is relatively good for the Earth, whose 
eccentricity is only about 0.0167. However, for Mars 
with an eccentricity in excess of 0.0937, the difference 
between the closest and farthest distance from the sun 
is almost 3/10 of an astronomical unit. 

The effects of the ellipticity of the Martian orbit 
were very apparent when the three-dimensional solu- 
tion for trajectory No. | was attempted. Since the 
difference of the mean longitudes of the two planets 
Earth and Mars is the same every 2.135 years, it would 
seem at first that departure dates that differ by this 
amount should produce essentially the same solutions. 
However, primarily because of the variation in the dis- 
tance from Mars to the sun, it was impossible to obtain 
a departure orbit for trajectory No. 1 with a departure 
velocity of 0.13 times the Earth’s orbital velocity before 
1964. For earlier dates, the point of contact with Mars 
occurred too near the aphelion of the planet. 

One final comment seems appropriate concerning the 
relative merits of the several trajectories considered. 
Although the departure velocities were all the same, 
the velocities with which the spaceship returns to 
Earth are seen to vary from 14,000 to well over 20,000 
it. per sec. These differences are of fundamental im- 
portance to the problem of re-entry of the vehicle into 
the Earth’s atmosphere. 


(3.4) A Class of Round-Trip Trajectories 


The particular solutions to the round-trip problem 
discussed in the previous section serve to illustrate 
some of the various types of possible reconnaissance 
trajectories. However, one should not infer that such 
solutions are isolated events. On the contrary, during 
appropriate seasons of the year, round-trip trajectory 
solutions exist as a continuous function of the departure 
time. It is of interest to determine how certain char- 
acteristics of these solutions vary with the time and 
velocity of departure. 

For this purpose, the particular class of trajectories, 
of which trajectory No. 7 is a representative member, 
was examined in full detail. Three departure velocities 
were selected—namely, 12,200, 13,000 and 13,800 ft. 
per sec. Then solutions corresponding to these veloci- 
ties were obtained at various instants of time, both 
earlier and later than December 4, 1964, which is the 
departure time for trajectory No. 7. To obtain a con- 
tinuous pattern of solutions, regular intervals of time 


24,000 


20,000} 


PASSING DISTANCE (MILES) 


oct NOV DEC Fee 
1964 
Fic. 3-10. Passing distance at Mars vs. departure date for 


various values of departure velocity. 


were selected and the time-of-flight estimates needed 
for the digital computer program were made sequen- 
tially, based on the time requirements of each previ- 
ously obtained solution. In this way, we were assured 
of producing only the trajectory solutions belonging to 
one particular class. The interval of departure times 
during which solutions were obtained using this pro- 
cedure extends from early in November, 1964, until 
late in March, 1965. Immediately prior to November, 
1964, and just after March, 1965, the configuration of 
the planets is unfavorable for a reconnaissance mission 
with the stated velocity specifications. 

In the accompanying figures, a few of the important 
elements of these solutions are plotted as a function of 
the departure time. The passing distance at the point 
of closest approach to the surface of Mars is shown in 
Fig. 3-10. During the first week in November, 1964, 
the point of departure from Earth and the point of 
arrival at Mars are nearly 180° apart. No solutions 
were obtained in this period because of the excessive 
velocity requirements out of the plane of the ecliptic. 
As time increases, the departure angle 6, decreases and 
round-trip solutions are possible. From late in No- 
vember to the early part of December. each of the con- 
stant-departure-velocity curves attains a minimum. 
Within ten days following the minimum, the required 
passing distance increases to such an extent that 
planetary contact becomes impossible. 
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Fic. 3-11. Velocity relative to Mars vs. departure date for var- 
ious values of departure velocity. 
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Fic. 3-12. Return velocity relative to Earth vs. departure date 
for various values of departure velocity. 


Early in January, 1965, solutions are again obtained; 
the required passing distance decreases steadily until 
sometime in the middle of February, when mathe- 
matical solutions are possible only by permitting the 
point of closest approach to fall below the Martian 
surface. During this period of time in February, the 
point of arrival at Mars and the point of return to Earth 
are nearly 180° apart. Later in February and early in 
March, the return angle 6g increases and physically 
realizable solutions are again possible. 

In Fig. 3-11, the velocity of the spaceship relative 
to Mars is plotted as a function of departure time. 
Similarly, the velocity relative to the Earth with which 
the spaceship returns to Earth is shown in Fig. 3-12. 
For the particular choice of departure velocities, the 
velocities relative to Mars vary from 18,000 to 27,000 
ft. per sec. during 1964 and from 25 000 to 29,000 ft. 
per sec. during 1965. The variations in the return ve- 
locity are much greater and, indeed, we observe the 
curious fact that the spaceship can return to Earth 
with a velocity smaller than the one with which it left. 

In 1964, all solutions require the spaceship to move 
through a heliocentric angular distance of less than 180° 
for the trip from Earth to Mars. During the return 
trip, the ship must orbit the sun once, traversing an 
angular distance greater than 540°. Roughly six 
months is consumed on the outbound leg, and two and a 
half years on the return voyage. In January and the 
first part of February of 1965, the outbound trip re- 
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Fic. 3-13. Angle between Martian and relative motion planes 
vs. departure date for various values of departure velocity. 
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quires a little more than a year and an angular dis- 
tance of some 250°. For the return trip, the ship must 
again orbit the sun and arrive back at Earth after 
approximately two years and an angular travel of 
less than 540°. Late in February and March the re- 
turn angular distance traveled is greater than 540°, 
while the outbound-voyage characteristics remain essen- 
tially unchanged. It is interesting to note that the 
total time of 3.2 years required for the trip does not 
vary by more than two months for any solution in the 
class. 

Concerning the manner of approaching Mars, the 
1964 solutions require the spaceship to cross the orbit 
of the planet heading away from the sun. The relative 
motion is such that the dark side of the planet is pre- 
sented to the ship during the greater portion of the 
period of contact. On the other hand, in 1965 the 
spaceship intersects the Martian orbit heading toward 
the Sun and the bright side of the planet dominates 
the field of view. In Fig. 3-13, a plot is shown of the 
angle between the normals of the Martian orbital plane 
and the plane in which the relative motion takes place. 
We see, for example, that in the first month and a half 
of 1965 the spaceship moves in a plane that is roughly 
perpendicular to the Martian orbital plane. 

In order to gain a clearer understanding of the char- 
acteristics of these round-trip reconnaissance trajecto- 
ries, a typical one has been mapped out in detail. The 
diagrams of Fig. 3-14 show various stages of the out- 
bound and return trip of a spaceship departing on 
January 25, 1965, with a velocity relative to the Earth 
of 13,800 ft. per sec. The orbits of the spaceship and 
of Mars are shown as solid lines when above the plane 
of the Earth’s orbit and as broken lines when they are 
below. 

Contact with Mars would be made on April 4, 1966, 
when the vehicle would pass within 4,880 miles of the 
planet’s surface. The orientation of the plane of rela- 
tive motion of the spaceship with respect to the Martian 
orbital plane is shown in Fig. 3-14c. In this illustra- 
tion, Mars is moving from right to left, the direction to 
the sun is indicated by the large arrow pointing toward 
the lower right, and the point of closest approach is 
shown by the small arrow. The primary effect of the 
Martian gravitational field is approximately to double 
the component of the spaceship velocity normal to the 
planet’s orbital plane. The result, as can be seen 
from the diagrams, is a rotation of the line of nodes of 
the spaceship orbit by some 60°. Furthermore, the 
inclination to the ecliptic of the orbital plane changes 
from 0.805° for the outbound path to 2.324° for the 
return path, and the period is reduced from 1.53 years 
to 1.49 years. Finally, the vehicle would return to the 
Earth on March 24, 1968, after a voyage lasting 3 
years and 58 days. 
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Design of Missile Bodies for Minimum Drag 


at Very High Speeds—Thickness Ratio, Lift, 


and Center of Pressure Given 


T. STRAND* 


Convair, A Division of General Dynamics Corporation 


Summary 


Newtonian flow theory has been used to develop a procedure 
for the design of minimum drag bodies of revolution having a 
given thickness ratio and center of pressure. 

It is shown that the optimum body shape is independent of 
lift. Center of pressure location, however, exerts a powerful 
influence on both the shape of the body and on the drag co- 
efficient at zero lift. 


Discussion 


t IS WELL KNOWN that when the speed of a missile 

becomes sufficiently high the shock wave lies very 
close to the body surface. In this case the local pres- 
sures can to a first approximation be calculated from the 
so-called Newtonian theory, which simply states that 
for the surfaces that ‘‘see’’ the flow the pressure co- 
efficient C, = 2 cos* (V, 2), where (V, m) indicates the 
angle between the free-stream velocity vector and the 
local normal to the body. Theoretically, the above re- 
sult is only valid for a gas at a free-stream Mach Num- 
ber of infinity and whose adiabatic exponent (y) is unity. 
Tests show that although the variation of the pressure 
along the body is accurately predicted by this formula, 
the level (the factor 2 above) is not correct. Recently 
a so-called ‘‘modified’’ Newtonian theory (that modifies 
the factor 2 above) has been suggested! and shows very 
good agreement with experimental results. According 
to this theory the stagnation point level is given by 


Comaz (Cy + 3)/(y + 1)] X 
{1 — [2/(y + 3)](1/M?)} (1) 


Tests show that this approach is valid from M = 2 or 3 
on up as long as the hypersonic similarity parameter (17 
times local slope) is somewhat larger than unity. 

In this discussion we shall assume that C,,,,, = 2 as 
given by the unmodified Newtonian theory. This does 
not invalidate any of the results obtained as the op- 
timum shapes are not a function of what constant is 
used for C,,,,,- It should be remembered, however, 
that all drag and lift values calculated here should be 
scaled down for any finite Mach Number by the ratio of 
Comar from Eq. (1)] to 2. 

Since all high-speed missiles have a blunt nose whose 
size is dictated by thermodynamic considerations, let us 
start by assuming that the nose is a hemisphere, and let 
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us first calculate the drag of the nose section. Here 
C, = 2 sin? @ = 2r’?/(1 + r’?) (2) 


where 6 is the local slope of the surface with reference to 
the x axis, r’ = dr/dx where r is the radius in a cylindri- 
cal coordinate system (Fig. 1), and x is the abscissa. 

If we let Ay be the radius of the hemisphere, the zero 
lift drag coefficient can be calculated from the following 
expression : 


0 
© (2 Ry?) f rr’ dx (3) 
— Ro 


Performing the integration, it is seen that Cp, = | 
which is in close agreement with the experimentally ob- 
served overall sphere value of 0.915 (Ref. 2), when the 
modified Newtonian correction is applied. (It can be 
shown that the skin-friction contribution in most cases 
is negligible for a sphere.) 

Next it is assumed that this nose drag can be added to 
the airframe drag, and that the hemisphere is small 
enough so that its contribution to the lift of the missile 
can be disregarded. 

We shall now restrict our attention to flight at small 
angles of attack. The angle of attack must be smaller 
than the smallest body slope in order to have the flow 
“see” all of the surface. The lift and drag can then be 
expressed as follows: 


L=N-—-C-a D=C+WN-a (4) 


where NV and C represent the normal and the chord 
forces, respectively. We shall now assume that the 
component of the chord force in the lift direction is of 
an order of magnitude smaller than the normal force. 
This puts an upper limit on the size of the nose radius 
as compared to the base radius. With this assumption: 
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Fic. 1. Coordinate system. 
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L@N, D=C+ N-a (5) 


In addition, it is assumed that the body slopes aft of 
the hemisphere are small compared to unity. The 
validity of this assumption must also be checked after 
the optimum body shape has been obtained. Then 
using Eq. (5) and the expressions for the normal and 
chord force coefficients developed in references 3 and 4, 
we have 


tage r’ dx => 2[1 (Ro |ag(rR,”) (6) 
0 


from which it is seen that C,, = 2/1 — (Ry Rp)*] inde- 
pendent of body shape. Ky, is the base radius. 

Let 7 be the pitching moment about the origin. 
Then 


M = f xrr’dx (7) 
0 
= 2[1 — (Vol./7Rg*) |ga(rRp*) (8) 
Here Vol. = (9) 
0 


It is seen that Cy, = 2[1 — (Vol./rR,?)|. These re- 
sults are identical to those obtained in reference 5 for 
inclined slender bodies of revolution using linearized 
theory. 

With the addition of the drag of the hemisphere 
found previously, the expression for the total pressure 
drag becomes 


1 
D = + [ + (3 2)L-a@ (10) 
0 


from which it is obvious that 
dCp dC,” = (3 4) [1 (Ro Rp)? | (1 1) 


for all bodies of revolution of the type under consider- 
ation. 

Next let us consider the problem of finding the body 
of revolution having minimum drag for a given thick- 
ness ratio, lift, and center of pressure. Problems of a 
similar type have been treated in reference 6. It is im- 
mediately obvious from the above results that the op- 
timum body shape is not a function of the lift. Hence, it 
is only necessary to minimize the integral in Eq. (10) 
while constraining the volume integral. This is then 
equivalent to specifying the center of pressure. Using 
the calculus of variation we form a new function /: 


[= [2 + dr? |dx (12) 
0 


where J is the constant Lagrangian multiplier. Let / be 
the integrand 


f=2rr'* + dr? (13) 
For a maximum or a minimum it is now necessary that / 
satisfy the Euler-Lagrange equation of the calculus of 
variation——namely, 


r'f,, — f = const. (14) 
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since the integrand functions are free of the independent 
variable. Substituting Eq. (13) into Eq. (14), we ob- 
tain: 

4r — dr? = C, (15) 


where C is a constant to be determined by the boundary 
condition. Defining K*? = A/C, and integrating Eq. 
(15) once more, it is easy to show that 


Kr 
x= (4 [o/(1 + p?)]' “dp (16) 
K is assumed to be positive. The negative and positive 
signs in the denominator of the integrand are used to 
take into account the fact that A/C, will change signs. 

When the negative sign is used, p is always less than 
unity, because we are interested only in the case of x in- 
creasing with increasing radius. 

When x = 1, r = Rg. Hence, from Eq. (16): 


= 1) [p/(1 + p?)]! “dog (17) 
KR, 


which, when substituted into Eq. (16), yields the equa- 
tion of the optimum body of revolutien: 


Kr KRB 
x= f [p/(1 + [p/(1 + p*)]' 
KR, KRo 


(18) 


The distance to the center of pressure is given by 


dCy, or Cra: 
= [1 — (Vol./#Rg?)]/[1 — (Ro/Ra)?] (19) 


Xe.p. 


The procedure is now for given values of Ry and KR, to 
choose a value of the constant AK and compute the 
body shape (Eq. 18) and the corresponding location of 
the center of pressure (Eq. 19). After a couple of trials 
the constant K for the body shape having the desired 
center of pressure location can be pinpointed. 

When K = 0 an indeterminancy results which cor- 
responds to the case of \ = 0; i.e., the case when no 
constraint except thickness ratio is put on the optimum 
body shape. It may easily be shown that in this case 


| Xop 
= + + 
) 1} + = + 4 — Pe 
: 
A) 2 a2 10 
+ 
4A 
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Fic. 3. Minimum drag of bodies of revolution with given 
thickness ratio and nose radius. (Note that drag coefficient 
does not include skin friction. ) 


= + (Ra? — (20) 


and that 
8/5) Ra) + (2/5)(Ro/ Ro)" 
[1 — (Ro/Rzs)*] 1 — (Ro/Rz)?] 


A few of the optimum body shapes for ky = 0.02 
and R, = 0.10 and their center of pressure locations 
are shown in Fig. 2. 

The pressure drag coefficient at zero lift as given by 
Eq. (10) is 


1 
Cp, = (Ro/ Rr)? + rr’ dx (22) 
0 


A graph of this coefficient as a function of the center 
of pressure location for the sample bodies of revolution 
is given in Fig. 3. 
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It is sometimes desirable to design bodies of revolu- 
tion which feature a minimum variation of center of 
pressure with lift (i.e., dCy/dC, ™ constant). In this 
case the center of pressure at low angles of attack should 
probably be located at the centroid of the area cut by 
a plane through the x axis. (The centroid of the area is 
the center of pressure location for a body of revolution 
in Newtonian flow at 90° angle of attack. Locating the 
centers of pressure at the same point at both a ~ 0° 
and a ~ 90° does not however automatically ensure 
no center of pressure travel for  interimediate 
angles of attack.) For the sample bodies of revolution 
it can be shown that when AK ~ 20 the center of pres- 
sure and the centroid are at the same location. The re- 
sulting body, which should have nearly constant center 
of pressure independent of lift, is almost identical to 
the third body from the top in Fig. 2. A drag penalty of 
ACp, = 0.001 is indicated for this body, however, com- 
pared to the absolute minimum drag body (K = 0). 
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Similar Temperature Boundary Layers 


DRAGUTIN STOJANOVIC* 
Unwersity of Beograd, Jugoslavia 


Summary 


Conditions for the existence of similar solutions are known for 
(a) two-dimensional, incompressible, steady and nonsteady 
laminar boundary layers and (b) three-dimensional, incom- 
pressible, steady, laminar boundary layers for a body of revolu- 
tion rotating in a fluid at rest or a body of revolution in a rotating 
fluid flow. Corresponding conditions for the existence of similar 
temperature boundary layers in both cases are given for constant 
and variable wall temperatures. The genera! conclusion is that, 
in all these cases, with or without viscous heating, and with con- 
stant wall temperature, conditions for the existence of similar 
velocity boundary layers are at the same time the conditions for 
the existence of similar temperature boundary layers. If the 
wall temperature is variable, the conditions for the existence of 
similar velocity boundary layers are at the same time the condi- 
tions for the existence of similar temperature boundary layers if 
the wall temperature varies as a power of the local free-stream 
velocity or surface velocity. Numerical solutions are given 
for the nondimensional temperature distributions function and 
the nondimensional temperature gradient at the wall for several 
Prandtl Numbers in the case of a rotating flow over an infinite 
plate at rest. 


Symbols 

é; = constants 

= specific heat 

f = factor depending on position for body of 
revolution 

F = nondimensional velocity distribution func- 
tion 

G = nondimensional velocity distribution func- 
tion 

g = gravity acceleration 

k = thermal conductivity 

L = characteristic length 

m,n = constants a 

p = pressure 

r = radius 

Re = Reynolds Number 

t = time 

r = temperature variable 

T @ = free-stream temperature 

i = wall temperature (constant cr character- 
istic) 

= coordinates 

u,v, W = velocity components in x, y, and =z direc- 
tions, respectively 

uy = local free-stream velocity for two-dimen- 
sional boundary-layer flow 

U = surface velocity of a body of revolution ro- 


tating in a fluid at rest or local free- 
stream velocity on a body of revolution 
in a rotating fluid flow 
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U. = characteristic velocity for a body of revolu- 
tion 

a, B, y, 6, Yo, € = constants 

= density 

v = kinematic viscosity 

= dynamic viscosity 

o = Prandtl Number 

én, ¢ = nondimensional coordinates 

¢, v = stream functions 

w = angular velocity 

0 = nondimensional temperature function 

) = function defined in text 


Similar Temperature Boundary Layers in Steady 
and Nonsteady Two-Dimensional Flows 


T IS KNOWN that in special cases the exact solutions 

for the two-dimensional boundary-layer equations! 
for steady, incompressible, laminar flow with variable 
free-stream velocity can be found (see Fig. 1). These 
“similar’’ boundary layers have the property that the 
velocity distributions in boundary layers differ in differ- 
ent positions only by one factor, which depends on 
position. Such similar solutions are also found for 
corresponding temperature distributions in boundary 
layers with constant and variable wall temperature. 
H. Schuh? has found the conditions for the existence 
of similar two-dimensional boundary layers in in- 
compressible, nonsteady, laminar flow with variable 
free-stream velocity. The related conditions for the 
existence of similar temperature boundary layers are 
found below. 
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Fic. 1. Coordinate system notations. 
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Starting from the equation for two-dimensional, 
incompressible, nonsteady, laminar boundary layers, 


(Ou Ot) + u(Ou Ox) + v(Ou/dy) = 
—(1/p) (Op/Ox) + v(07u/dy?) (1) 


the equation of continuity, 
(Ou/Ox) + (Ov/Oy) = 0 (2) 
the conditions at the boundaries, 
u=v=0, (x,t), yoo (3) 
and the free-stream velocity-pressure relationship, 
—(1 p)(Op Ox) = (Om,/Ot) + (4) 
Schuh indicated that with solutions of the form 
u = u,:F'(n), n = y/N(x, t) (5) 
Egs. (1), (2), and (4) are transformed into 


— (N*/v) (Om/O0x) (F”? — F’F) + 
(N v) (ON /Ot)F’n + (uN /v) (ON/Ox)F"F — 
(N? (Om, Ot) F’ + (N?/v) (Om,/Ox) + 
(N?/ vis) (Om; /Ot) =() (6) 


From Eq. (6) we see that similar solutions exist if the 
relations 


(N?/vu1) = (7) 
(N?/v) (0u;/Ox) = Ce (8) 
(N/v) (ON/Ot) = 3 (9) 
(u:N/v) (ON/Ox) = (10) 


are constants ¢, C2, C3, Cs. Then the partial differential 
equation (6) becomes the ordinary differential equation 
(11), 
+ + — — — 
aF’ (11) 
with the conditions at the boundaries, 


To find the corresponding conditions for similar 
solutions for thermal boundary layers, we start from 
the energy equation for incompressible two-dimensional 
boundary layers. If the fluid properties are constant 
this equation is 


pgc,[(OT/Ot) + u(OT/Ox) + v(OT/dy)] = 
k(0?T /Oy?) + w(Ou/dy)? (13) 


We can find the general solution in the form 
T = d(x, t)-6(n)-Tw, To = 0 (14) 
with the boundary conditions 


T = 0,(x,1)-To, y= 0, 0=1, (15) 
T— 0, y>o,6>0, of 


It is clear also that we can find by superposition solu- 
tions of the form 


T = A(x, t)-S(n) + T,8 (16) 
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with the boundary conditions 


S=0, 6=1, 
S=0, 06=0, 


n=0 | 


(17) 


but, because this form of the solutions adds nothing 
special to the generality of the conditions for existence 
which we are looking for, we use only the form (14). 

With the transformations for velocities, with condi- 
tions on the velocity boundary layers, Eqs. (8), (9), 
(10), and, with transformations, Eq. (14), the energy 
equation (13) can be transformed to 


A” + a(Cy + C,) Fo’ + onc;0" 
(pgc,/k) (N?/d:) — 
(pgco/k) + 
(u/RT (17/01) = 0 (18) 


With (pgc,/k) = C5 (19) 
(pgc,/k) (N? (Ov; Ot) = (20) 
(u/Tk) = (21) 


the energy equation can be written in the form 


6” + a(cs + cs) FO’ + — — 
+ oF" = 0 (22 


From this equation it can be seen that in similar 
boundary layers there can also exist similar temper- 
ature boundary layers if cs, cs, and c; are constants (the 
case in which c;, cs, and c; are functions is not con- 
sidered). In that case, Eq. (20) is an ordinary differ- 
ential equation and the conditions (19), (20), and (21) 
are the general conditions for the existence of similar 
temperature boundary layers. The form of this 
equation suggests the influence of similarity factors of 
the velocity boundary layers, 4, C2, cs, cs, and of the 
thermal boundary layers, C5, Cs, C7. 

Now we are going to observe three special cases of 
these conditions. The first special case is that in which 
wall temperature is unvariable and viscous heating is 
neglected. In this case, 


and Eq. (22) is simplified to 
AY + + cs) FO’ + = 0 (23) 


As this equation has only constants ¢2, ¢3, and ¢, 
which are also the conditions for the existence of 
similar velocity boundary layers, it is clear that in this 
case similar temperature boundary layers always exist 
and the conditions for the existence of similar velocity 
boundary layers are at the same time the conditions 
for the existence of similar temperature boundary 
layers. 

The second special case is that in which wall tem- 
perature is variable and viscous heating negligible. 
In this case c; = 0 and conditions for the existence (19) 
and (20) can be written in the more convenient form 


(1/2) (1/t1) (Ou,/Ox) = (1/81) (081 /Ox) (24) 
(o/Cs) (1/01) (25) 


(1/e1) (1/11) 
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The general solution of these equations is 
3, = au," (26) 
and in this case we can find 


Cs = ONCo, Ce = Onc, (27) 


Then Eq. (22) is transformed to 


+ + cx) Fo’ + 
onc,F’6 — = 0 (28) 


and from this equation it can be seen that if the wall 
temperature varies as a power of the local free-stream 
velocity the conditions for the existence of similar 
velocity boundary layers are at the same time the 
conditions for the existence of similar temperature 
boundary layers. It is interesting to note that, in the 
case of variable wall temperature, in transformed Eq. 
(28) there also exists the similarity constant, c, of the 
velocity boundary layers. 

In the third case, if we consider the complete Eq. 
(22), the similarity conditions for thermal boundary 
layers, Eqs. (19), (20), and (21) give 


and the complete transformed equation for thermal 
boundary layers is 


20¢,0 — moF”"? = 0 (30) 


We can say that with viscous heating the conditions for 
the existence of similar velocity boundary layers are at 
the same time the similarity conditions for the thermal 
boundary layers. 


Similar Temperature Boundary Layers on a 
Body of Revolution in a Rotating Fluid or on a 
Rotating Body of Revolution in a Fluid at Rest 


It is known that for some bodies of revolution similar 
solutions of the boundary-layer equations can be found 
if the flow is steady, axially symmetrical, incompres- 
sible, and laminar (see Fig. 2). For this case, starting 
from the boundary-layer equations of W. D. Hayes, 
Th. Geis* has found that boundary-layer equations 
for a body of revolution in a rotating fluid flow are 


vu, + wu, + uv(r’/r) = vu, (31) 
w, + we, — urr’/r) = —U%(r'/r) + wz. (32) 
Vy + wr + vo(r’/r) = 0 (33) 
with conditions at the boundaries 
u=v=w=0, 2=0 
u—>U, 


and that boundary-layer equations for a body of revolu- 
tion rotating in a fluid at rest are 


+ wu, + ur(r’/r) = vuz, (34) 


Wy u*(r’ ‘r) = We (35) 


Z 


axis of revolution 
x 


Fic. 2. Coordinate system notations. 


v% + w, + v(r’/r) = 0 (36) 
with the boundary conditions 


or = U, v=w=0, s=0 


u—~>O0, 27> & 
Geis has transformed these equations with 
W=¢, W= (37) 
t=x/L, n=y/L, ¢ = 2VRe/L-f(y), 
Re = U.L/v_ (38) 
F(n, §) = 8) [WRe/L- U(y)-f9)| \ (30 
G(n, = 2) [WRe/L- 
to (a) for a body of revolution in a rotating fluid, 


+ oGF" + BF'G' = 0 
+ aGG" + + - 1) =of 


with the boundary conditions 


F= F'’=G=G' =0, 


i= 


and (b) for a body of revolution rotating in a fluid at 
rest, 


+ aGF" + BF’G’ = 0 (42) 
+ aGG" + 7G" + 6F”? = 0! 
with the boundary conditions 
F=G=G' = 0, 
F’—+0, G’>0, (> (43) 


Conditions for the existence of similar velocity 
boundary layers are that the coefficients, 


a= —B+(1/2)0(f)y 8 = 
r = (r/L) 


a, 8, y, and 6 are constants. 

To find corresponding conditions for the existence of 
similar temperature boundary layers for these cases of 
fluid flow we start from the energy equation. If 
physical properties of fluid are constant and all bound- 
ary-layer simplifications are made, this equation can 
be written in the following form: 


pgc.|v(OT /Oy) + w(OT/dz)] = + 
u[(Ou/Oz)? + (Ov/dz)*] (45) 
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FIG. 3. 


constant wall temperature. 


with the boundary conditions 


T = 

With 
T = Ty (47) 


and with all similarity conditions from the velocity 
boundary layers, the energy equation can be trans- 
formed to 


6” + o(a +B +6 — y)GO’ — od(r/r,) X 
+ (u/k) (1/T») x 


+ G") = 0 (48) 


with the boundary conditions 


6=1, ¢=0, 00, (49) 


Conditions for the existence of the similar temper- 
ature boundary layers are that, as well as a, 8, y, 6, 
being constants, 


(r/r_) = 1(O/U,) = € (50) 
(u/k) (U.2/Tw) (07/01) = (51) 

are constants also. 
Again we consider some special cases. The first 


special case is that of an invariable wall temperature 
Then the constants « = ¥y = 0 and Eq. (48) is simpli- 
fied to 


6” + ola +B +6 — =0 (52) 


The conditions for the existence of similar velocity 
boundary layers are at the same time the conditions for 
the existence of similar temperature boundary layers. 
If constants have the values 

it 

+t 


a= -2, 6= 
B=+2, y= 
then this case is that of an infinite plate rotating in a 
fluid at rest or an infinite plate at rest in a rotating fluid. 
The corresponding exact solutions for temperature dis- 
tributions were found earlier for the case of the rotating 
plate. Here we give the corresponding solution for an- 
other interesting case of a plate at rest in a rotating 


(53) 
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Nondimensional temperature distribution for rotating 
flow over an infinite plate at rest for various Prandtl Number and 
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fluid. With numerical values for velocity functions 
from reference 5, the numerical values for temperature- 
distribution functions are computed, since in this case 
the equation is simplified to 


6” — 2oG0’ = 0 (54) 


and the solution can be found from 


er £ Gdt a) (55) 
0 0 


for various Prandtl Numbers. Figs. 3 and 4 give these 
values for the temperature-distribution function and 
the temperature gradient at the wall. 

The second special case is that of variable wall tem- 
perature and negligible viscous heating. From the 
similarity conditions (50) we can find that, for the 
existence of similar temperature boundary layers in 
similar velocity boundary layers the wall temperature 
must vary as 


3d, = AU" (56) 
then Eq. (48) is transformed to 
0” + o(a +8 +6 — — aynG'6 = 0 (57) 


The third special case is that in which viscous heat- 
ing exists and the wall temperature is invariable. 
From the conditions of existence (50) and (51) we can 
find 

= (U2/T,mgc,) (58) 
and Eq. (48) is transformed into the ordinary differ- 
ential equation 


6" + o(a + B+ 6 — — + 
ma(f’? + = 0 


(Continued on page 591) 


(59) 


Fic. 4. Temperature gradient at the wall for rotating flow 


over an infinite plate at rest for various Prandtl Number and con- 
stant wall temperature. 
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Similitude of Hypersonic Real-Gas Flows Over 
Slender Bodies With Blunted Noses’ 


HSIEN K. CHENG* 


Cornell Aeronautical Laboratory, Inc. 


Summary 


On the basis of the hypersonic small-perturbation theory, the 
laws of similitude for hypersonic inviscid flow fields over thin or 
slender bodies are examined, and the restrictions to ideal gases 
with constant specific heats and to bodies with pointed noses are 
removed. Only steady plane or axisymmetric flows are con- 
sidered. 

Inspection of the governing system of equations shows that a 
similitude law exists for flow fields, under local thermal equilib- 
rium, having the same free-stream atmosphere. For flows of 
ideal gas with constant specific heats, the requirement of the same 
free-stream atmosphere—i.e., the same composition, pressure, and 
density—can be replaced by the requirement of the same ratio of 
specific heats, 

For flows over blunted wedges or cones, special laws of simili- 
tude can be obtained. 

Application of the similarity rules is examined for the case of 
hypersonic flows of an ideal gas with y = 1.40 over flat plates 
with blunt leading edges, and for the case of equilibrium air 
flows over wedges. The possibility of simulating nonequilibrium 
flows over slender or thin bodies is also pointed out. 


Symbols 


x,y = linear normal coordinates for plane and 
axisymmetric fields; x gives the distance 
from the origin along the axis parallel to 
the uniform free stream, and y the dis- 
tance from the axis 

v = a number, vy = 0 and 1 for plane and axi- 
symmetric flows, respectively 


Ue the undisturbed free-stream velocity 

u,v = the perturbation velocity components in 
the x- and y directions, respectively 

Gy, Cy = specific heats of the gas at constant pres- 
sure and at constant volume, respec- 
tively 

Y = C,/C, the ratio of the specific heats 

d.0;7T = the pressure, density, and temperature of 


the gas, respectively 
h, e, s = the specific enthalpy, the specific internal 
energy, and the specific entropy of the 
gas, respectively 
the concentration (either in mass fraction 
or in number fraction) of the species 7 
in a gas mixture 
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Subscripts 


@+(1/M.) = 


(pol a/pao)(p/pT) = the ‘compressibility 
factor’’ of air 

the equilibrium speed of sound in the free 
stream 

the free-stream Mach Number 

the lateral extent of the hypersonic flow 
field and the coordinates of the bow 
shock wave, respectively 

the thickness of the layer of flow which has 
passed through the blunt-nose region, 
see Fig. 1 and Appendix 

a_parameter characteristic of the flow de- 
flection angles 

a typical dimension of the body in the 
streamwise direction, see Eqs. (10) and 
(19) 
a parameter used in the 
hypersonic small-perturbation theory 
the function which defines the body sur- 
face, see Eq. (5) 

a dimensionless parameter, see Eqs. (10) 
and (19) 

the angle of attack of the wedge surface, or 
the half cone angle 

variables and functions defined by the 
similitude transformation, Eqs. (10) 
and (11) 

the nose drag 

the nose drag coefficient defined according 
to Eq. (17) 


= Lift/(1/2)p.U°L, lift coefficient 


Pitching Moment/ [(1/2)p. U?L?], moment 
coefficient 

Drag Force/|(1/2)p.U*%Area)], see Eqs. 
(22) and (23) and the explanation which 


follows. Note that the nose drag has 
been included 

be — po = pressure increment on the 
surface 


function of the variables &, .. ., ete., 
involved 

a typical dimension of the blunt-nose re- 
gion in the x direction, see Fig. 1 

a typical dimension of the thickness of the 
nose in the y direction, see Fig. 1 and 
Eq. (17) 

the stream functions, see Fig. 1 and Ap- 
pendix 


© pertaining to the free-stream condition 

N pertaining to the region of the blunt nose 

* pertaining to the thin layer of flow which passes through. 
the blunt-nose region 

B pertaining to the flow over the body surface 

0 pertaining to the sea level atmospheric condition 


Introduction 


ys SIMILARITY LAWS of hypersonic inviscid flows 
of an ideal gas with constant specific heats over 
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thin bodies were first given by Tsien.’ In Tsien’s orig- 
inal small perturbation theory for hypersonic flows, 
only the class of irrotational, plane, or axisymmetric 
steady flows was considered. Subsequently, Hayes’ 
pointed out the equivalence between the field of a hyper- 
sonic flow as observed from the thin or slender body and 
the transient flow field in a transverse plane; and the 
laws of similitude may be regarded as being character- 
istic of the transient cross-flow problem. The restric- 
tion to the class of plane or axisymmetric irrotational 
flows implicit in the earlier work of Tsien is thus seen 
to be unnecessary. 

Further studies of the law of similitude have been 
made by many authors (for example, see references 
3-10) on the basis of the governing partial differential 
equations and the boundary conditions for flows of 
ideal gases with constant specific heats. These in- 
vestigations confirm that similitude applies only among 
flow fields having the same ratio of specific heats, y, 
as well as the same value of the hypersonic parameter 
Mr. 

For flight at hypersonic speed in the earth's atmos- 
phere, because of the high temperature level encoun- 
tered, the real gas effects associated with internal exci- 
tations, such as vibration, dissociation, and ionization, 
play an important role in determining the gas proper- 
ties in the field.'~*! In particular, the assumptions 
that C,, C,, and p/pT are constant, which are implicit 
in the previous studies, appear to be inadequate. A 
question then arises as to the specific conditions under 
which the simulation of hypersonic flows of real gases 
is possible. Moreover, the similitude of Tsien and 
Hayes is not applicable to a class of flow fields of con- 
siderable practical interest—namely, the flows over 
thin or slender bodies which possess small blunt leading 
edges or noses. 

In the present inquiry, the governing differential 
equations and boundary conditions for the inviscid 
hypersonic flows over thin or slender bodies are used to 
extend the laws of similitude to real gases and to bodies 
having blunt noses. 

In the following, only the similitude of steady flow 
fields in which local thermal equilibrium prevails will be 
examined. The possibility of simulating real-gas 
flows with relaxation effects will also be pointed out. 
In order to derive the similitude laws in the most gen- 
eral form, the physical nature of the gas will be unspeci- 
fied, except that it is a nonviscous continuum. 


Basic Assumptions 


The assumptions underlying the following inquiry 
are: 

(a) The flow velocity is much greater than the speed 
of sound in the free-stream atmosphere. 

(b) The local flow deflection angle (in radians), as 
measured from the undisturbed flow, is much smaller 
than unity except perhaps in the nose region. 

(c) The free stream is uniform in velocity and 
homogeneous in the gas properties. — 
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(d) The excitation levels of all degrees of freedom of 
the gas molecules are in thermal equilibrium with the 
local pressure and density (or temperature). 

(e) The flow process is particle-isentropic, except 
when crossing the shock-wave surface. 

Assumptions (a) and (b) can be written as 


> 1, (hypersonic) 
6<1, (small perturbation) 


The parameter 6 controls and characterizes the magni- 
tude of the flow deflection angles measured from the 
free-stream direction; viz., 


= 


Hence, the phrase “‘small perturbation’’ used here refers 
only to the velocity field, and no restriction is to be 
placed on the magnitude of the variations in pressure 
and in density. In this paper, the effect of the blunting 
of the noses of thin or slender bodies will also be con- 
sidered. In such cases, the small-perturbation assump- 
tion is certainly not valid in the nose region, and the 
analysis will be applied only to regions sufficiently far 
downstream of the nose. 

Under thermal equilibrium, the thermodynamic prop- 
erties of a reacting gas mixture characterizing the state 
are determined by any pair of independent variables, 
such as p and p, provided the initial composition at a 
standard state is known. These properties have been 
tabulated for air and can be found, for example, in 
references 11-13. Accordingly, one can write, under 
assumption (d), the temperature, specific entropy, 
specific enthalpy, as well as the concentration of each 
species in the gas mixture, in the form: 


etc. 


The standard state is taken as that of the free-stream 
atmosphere described by p., pa, and X,.. The 
symbol f is used to represent a function of the variables 
involved, different for each of the quantities 7, s, ete. 
The X;, refer to concentrations of atoms and molecules 
of various species as well as the ions and the electrons. 

The particle-isentropic assumption (e) implies that 
the analysis can only be applied to the ‘‘nonviscous 
flow’’ regions where the transport properties of the gas 
are not essential in determining the fields of velocity, 
pressure, and density. 

The following section is concerned with the orders 
of magnitude of the streamwise component of the per- 
turbation velocity and the lateral extent of the per- 
turbed flow field, and with the pertinent differential 
equations and boundary conditions, restricting con- 
sideration to the pertinent flow field sufficiently far 
downstream of any blunt nose. In order to simplify 
application (in a subsequent section) to bodies having 
blunt noses, the present study is confined to plane and 
axisymmetric steady flows. An analysis similar to the 
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following can, in fact, be carried out for the general 
three-dimensional case.** 


The Governing Differential Equations and the 
Boundary Conditions 


A Formulation of the Small-Perturbation Theory of 


Hypersonic Flows 


Examination of the Rankine-Hugoniot relations for 
an oblique shock wave deflecting an hypersonic stream 
by a small angle reveals that the following orders of 
magnitudes hold behind the shock: 


v/U, = 0[6| 
= 0[6+ (1/M.)] = Ole], (1) 
u/U. = Ofe?| 


In the following, these estimates will be applied to the 
entire flow field between the bow shock wave and the 
body surface, excluding the region at the vicinity of a 
blunt nose. The justification will be seen from the 
consistency among the supposition Eq. (1), and the 
resulting system of equations governing v/ L’.., 6 L, and 
u/U.. as well as the basic assumptions (a)—(e). 

Assuming now the orders of magnitudes given in Eq. 
(1), the system of differential equations governing the 
motion of a particle-isentropic, plane or axisymmetric 
steady flow is reducible to 


[U..(0/Ox) + v(0/dy) |v = —(1/p) (Op/dy) | 
[U..(0/O0x) + v(0/dy) |s = 0 (2) 
[U..(0/Ox) + v(0/Sy)]o = —p[(dv/dy) + 


[U.(0/0x) + v (0/dy) = —(1/6) (Op/Ox) (3) 


where the symbol v is a number which is to be taken as 
zero for plane flows and as unity for axisymmetric 
flows. The specific entropy s, because of the thermal 
equilibrium assumption (d), is a known function of p 
and p. Thus the three equations (2), together with 
appropriate boundary conditions, are sufficient in 
principle to determine the fields of p, p, and v. The 
foregoing equations incorporate the order of magnitude 
estimates of Eq. (1); in particular, the terms associated 
with the differential operator u(0/0x) as well as the 
terms p(Ou/Ox) have been dropped, whereas the non- 
linear terms associated with the differential operator 
v(0/Oy) are retained. 

The streamwise perturbation velocity « characteriz- 
ing the three-dimensional effect may be determined 
from Eq. (3) after the equation system for p, p, and 7 
has been solved. In fact, Eqs. (2) and (3), together 
with the uniform free-stream condition (c), yield the 
Bernoulli relation from which u can be directly calcu- 
lated in terms of p, p, and v: 


2(u/U.) = —(v/U.)? — [(a—h.)/(1/2)U (4) 


where the specific enthalpy h, under thermal equilib- 
rium, is a given function of p and p. 

To the same degree of approximation, the boundary 
condition on the body surface y = Y(x) is 


g9= U.Y'(x) (5) 


and the Rankine-Hugoniot relations across the shock 
front vy = 6(x) are 


piv — = —p.U,6’ 
pb + — U.6’)? = p.U.*(8')? + p (6) 
(1/2) — +h = (1/2)U 


where 6’ = dé /dx and the shock surface y = 6(x%) may 
not be known a priori. In arriving at Eqs. (6), the re- 
lations Y’(x) < 1 and 6’(x) < 1 have been used, which 
are consistent with the small perturbation assump- 
tion. 

The above system of equations, Eqs. (2)-—((j), to- 
gether with the thermodynamic relations giving s and 
h in terms of p and p, will be sufficient under assump- 
tions (a)—(e), for determining the flow field. Among 
the exceptional cases are the flows over thin or slender 
bodies with blunt noses, of which a special study will 
be made in a subsequent section. Eqs. (2)—(6), except 
for the use of the unrestricted forms for the entropy and 
the enthalpy, are essentially the same as those used 
previously by Goldworthy® " and Van Dyke® in their 
formulations of the hypersonic small-perturbation 
theories for the ideal gases with constant specific heats. 

These equations have been reduced from the exact 
equations of motion of inviscid flows through the use 
of the orders of magnitude estimates given in Eq. (1). 
The error resulting from this simplification in each of 
the equations is generally of an order which is e* higher 
than one of the remaining terms in the equation. Con- 
sequently, the corresponding error in the solution is of an 
order ¢? higher than the solution itself. In this con- 
nection, it should also be pointed out that, inasmuch as 
the small-perturbation assumption in the velocity field 
can be satisfied, the nonlinear hypersonic small-per- 
turbation theory and the usual linear theory of Prandtl 
and Glauert are not mutually exclusive. Examination 
of the governing equations given above and the solu- 
tions of the linear equations shows that the domain ot 
validity common to both theories occurs for \J.. > 1, 

Before passing, the equivalence of the problem pres- 
ently considered and the transient problem in the cross- 
flow plane, as pointed out by Hayes,? may be noted. 
The equations governing the fields of p, p, and v which 
one would observe from a frame of reference at rest 
with respect to the undisturbed fluid can be obtained 
from Eqs. (2), (5), and (6) by substituting ¢ for x U.. 
Since the perturbation velocity component in the 
streamwise direction « does not appear in these simpli- 
fied equations, the transformed system of equations 
can be identified as the system which governs the un- 
steady motion of the same gas in a space of one less 
dimension. Therefore, within the small-perturbation 
approximation, the problem of an hypersonic inviscid 
gas flow over a thin or slender body, as far as p, p, and 
v are concerned, is identically the same as the problem 
of the transient motion of the same gas in a cross-flow 
plane produced by the unsteady movement of the body 
surface as it passes through this plane. 
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Pep 
(oc). EXAMPLE OF BLUNT-NOSED THIN OR SLENDER 
BODIES FOR WHICH THE DOWNSTREAM INFLUENCE 
OF THE NOSES MAY BE IGNORED. 


(b). EXAMPLE OF BLUNT-NOSED THIN OR SLENDER 
BODIES FOR WHICH THE DOWNSTREAM INFLUENCE 
OF THE NOSES MAY NOT BE IGNORED. 


>? Zao 

pow sHock WAVE 

BODY SURFACE, =O 


— 


(c). DESCRIPTION OF FLOW REGIONS AROUND A 
BLUNT-NOSED THIN OR SLENDER BODIES. 


Illustration of hypersonic inviscid flew fields around a 
blunt-nosed thin or slender body. 


Fic. 1. 


Application of the Hypersonic Small-Perturbation Theory 
to Flows Over Blunt-Nosed Bodies 


The small perturbation assumption 6 < 1 assumed 
in the foregoing analysis does not hold in the vicinity 
of the nose of a blunt-nosed body. In principle, the 
flow field sufficiently far downstream of the nose can 
still be treated by the small-perturbation theory—.e., 
by Eq. (2)—(€)—provided that the body is sufficiently 
thin or slender. In the literature,**~* examples of 
blunt-nosed bodies can be found for which the small- 
perturbation theory of hypersonic flows can be readily 
used to determine the solution for the flow field down- 
stream of the nose. Though not completely general,* 
the class of nose shapes represented by these examples 
appears to be general enough to cover a range of con- 
siderable practical interest. The blunt-nosed bodies 
referred to hereafter in this paper will be restricted to 
this class for which small-perturbation theory applies 
downstream of the nose. Typical examples are illus- 
trated in Fig. 1. 

For the purpose of examining the downstream in- 
fluence of a small but blunt nose, one may consider the 
flux of momentum and energy across an appropriately 
chosen control surface enclosing a portion of the body 
and arrive at an integral relation: 


5 (x) 
Dix) = [ee + (1/2)¢0"] (7) 
x 
where D(x) is the total inviscid drag of the body up toa 
* For example, this class of blunt noses does not include the 


type V(x) ~ x° where o < 2/(3 + v) for which the small-pertur- 
bation theory does not readily yield a solution.*4 
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distance x from the nose. In arriving at the above 
integral relation, the small-perturbation assumption jis 
applied only at the station x. Eq. (7), therefore, may 
be used when there is a blunt nose somewhere upstream 
of station x, so long as the small-perturbation assump- 
tion is valid, at least at x. Thus, downstream of x, 
the solution can presumably be carried forward by using 
the small-perturbation equations previously developed. 
Assuming that the dimensions of the blunt nose are 
small, one may apply the small-perturbation theory 
downstream of a point which is so close to the nose 
that Eq. (7) takes on the limiting form for x > 0: 


5(x) 
Dy = lim ple + (1/2)v?] (2ry)" dy 
J ¥(x) 
where Dy represents the nose drag. Eq. (8) will be 
used, together with Eqs. (5) and (6), as a boundary 
(initial) condition to determine the solution of the 
differential equations Eq. (2). 
In order to see more clearly the extent to which the 
blunt nose may affect the flow field downstream, one 
may write 


D(x) = Dy + bsa(dY/dx)dx (9) 


The range of integration with respect to x in the above 
pressure-drag integral does not include the “nose re- 
gion.”’ It is evident, from the above equation and 
Eq. (7), that the nose effect can be disregarded if the 
body is such that the pressure-drag integral beconies 
much larger than Dy. To determine solutions for this 
type of body, Eqs. (2)—(6), together with the knowl- 
edge of the thermodynamic properties of the gas, will 
be sufficient. Examples of blunt bodies for which the 
nose drag may be neglected are the slender, or thin, 
‘power-law bodies’ treated by Lees and Kubota,” 
the effective body resulting from the displacement 
effect of an hypersonic viscous boundary layer treated 
earlier by Stewartson,”* and blunted wedges or cones 
when the wedge or cone angles are sufficiently large. 
Fig. 1(a) illustrates bodies of this type. On the other 
hand, if the nose drag Dy contributes considerably to 
the total D(x) it is apparent from Eqs. (7) and (9) that 
its influence will be quite large. To determine the 
downstream flow field, a boundary (initial) condition, 
Eq. (8), has to be used in addition to the boundary 
conditions on the shock and on the surface. Examples 
in this class are the ‘‘flat-plate afterbodies’’ and the 
“cylindrical afterbodies"’ treated by Lees and Kubota”! 
and Cheng and Pallone.”* 

In terms of the transient cross-flow concept of Hayes, 
Eq. (7) may be interpreted simply as an overall energy 
balance of the cross-flow field at a certain time ¢ = x/U). 
The right-hand side of Eq. (7) gives the total energy in- 
crease of the field, and the left-hand side gives the total 


} The term poe has been dropped in the integrand. The rea- 
son is seen from the fact that in general pe/p.oeo = O[p/po], and 
that for flows over blunt-nosed bodies, 6’ > ©,asx—~0O. Hence, 
from the Rankine-Hugoniot relations, pe/pseo ~ p/po ~ 
M ~*(5')?—> ~,asx—0. 
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energy delivered to the cross-flow plane. The effect 
of a small but blunt nose on the transient cross-flow 
field at time ¢, in view of Eq. (9), is therefore equivalent 
to that of a blast wave releasing an amount of energy 
Dy at the instant ¢ = x/U. = 0. In fact, Eg. (8) 
specifies such an initial condition. 
_ The cross-flow energy relation, Eq. (7), has been ap- 
plied by Lees and Kubota*‘ to determine the solutions 
for the flow fields over the flat plate and cylindrical 
afterbodies. The same relation has been used earlier 
by Taylor,* Lin,” and Sakurai® for the mathemati- 
cally equivalent problems of blast waves. 

As pointed out previously, the “‘initial condition” 
Eq. (8), rather than the integral relation Eq. (7), will be 
used in the present formulation of the problem. As 
far as determining the solution for the system of equa- 
tions, Eqs. (2)—(6), is concerned, Eq. (7) and Eq. (8) 
are equivalent. For Eq. (7) in essence is an integral 
of the differential equations, Eqs. (2), which is con- 
sistent with the boundary conditions on the shock and 
on the body surface. The only additional information 
embodied in Eq. (7) relates to the energy release Dy 
to the transient cross-flow field at 4 = x/U., = 0 which 
is just as well specified by Eq. (8). Using Eq. (8), 
together with Eqs. (5) and (6), as boundary conditions 
isin accord with usual practice in the specification of 
boundary value problems. 

While the magnitudes of the errors implicit in the 
hypersonic small-perturbation theory have been pointed 
out previously, it should be noted that, in cases in- 
volving blunt noses, the magnitude of the error in the 
solution for the density or the temperature increases 
without limit approaching the body surface, giving a 
vanishing density or infinite temperature there. * 
The failure of the theory to provide a meaningful 
answer for the density or temperature distribution near 
the surface of a blunt-nosed body is related to the fail- 
ure of Eqs. (2)—(6) to deal with the flow field in the 
blunt nose region where the flow deflection angles can 
no longer be taken as small. The error (in entropy) 
resulting from violation of the small-perturbation as- 
sumption in the nose region persists in the flow down- 
stream of the nose region. The region of (uniform) 
validity of these solutions therefore excludes both the 
neighborhoods of the blunt nose and that of the body 
surface. The lateral extent of the singular region 
near the body surface 6 is, however, small in com- 
parison with the lateral extent of the hypersonic 
field 5, provided the small-perturbation assumption can 
be locally satisfied. In this case, the pressure vari- 
ation across the layer is negligible, and the surface pres- 
sure predicted by the small-perturbation theory is 
sufficiently accurate. A more detailed discussion of 
these points is given in the Appendix. 

The hypersonic small-perturbation theory, by the 
use of Eqs. (2)-(6), and (8), and the thermodynamic 
relations giving s and / in terms of / and p, has been 
formally generalized to include the effects of real 
gases as well as blunting of the noses. On the basis of 
the above system of equations, the similitude of steady 
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hypersonic flows will be examined in the following 
section. 


The Similitude of Hypersonic Flow Fields 


General Similitude 


For the purpose of studying similitude, an affine 
transformation is introduced 


v= p=pob, p 
where L is a reference length and 7 an arbitrary param- 
eter which will be seen to be a thickness parameter of 
the body. Under the thermal equilibrium assumption 
(d) we can also express the thermodynamic variables 
s, h, ete., in the form: 


(11) 


With these new variables, the system of differential 
equations, Eqs. (2), is transformed to 


ja = —(0p/d5) (12) 

[(0/Ok) + = + v(8/5) 

+ 0(0/0%) = 0 

The boundary condition on the body surface, Eq. (5), 

becomes 


= 6 = (13) 
The condition across the shock front, Eq. (6), becomes 


= 4(2), 
pls — = —6'(2), 
p (a..*p./p a) (M ..r)*plé |? 
(a (Mr)? [6 — +h = 
(a .2/2h..) (M or)? +1 
The initial energy condition, Eq. (8), may be written 


=lim [(a..2/2e.) X 


Y (2%) 


(M + pe] (15) 


The streamwise component of the perturbation velocity 
becomes 


2a = — (h./2a.”) — 1)/(M.7)*] (16) 
Eqs. (11)-(16), in principle, determine the solution 


for p, p, 5, h, 3, and 6 in the variables Z and 7. The 
following parameters appear in the system: 


Mat, 


The group of thermodynamic parameters pertaining 
to the ambient condition can simply be replaced by 
Po, po and X;,., which specifies the free-stream 
atmosphere under thermal equilibrium. The second 
parameter may be reduced to 
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where k is the nose drag coefficieat defined as 
k = Dy/(1/2)p «*(d/2)[w(d/2) (17) 


Hence, the solution to Eqs. (11)—(16) is completely de- 
termined by the following parameters: 


This solution, according to the affine transformation, 
Eq. (10), yields a family of ‘‘similar’’ hypersonic flow 
fields for various combinations of LU’ .., 7, k, d, and L. 

A law of similitude for hypersonic flows is now ob- 
tained: For a steady hypersonic flow field satisfying 
the assumptions (a)-(e), there exists a family of 
“similar’’ flow fields corresponding to the same values 
of M.rand (d/L) and the same 
free-stream atmosphere. Furthermore, the body shapes 
(away from the small blunt nose) are described, in the 
dimensionless coordinate system, by the same equa- 
tion: 


= P(x/L) (19) 


The similar flow fields are given by the same solution to 
Eqs. (11)—(16) in the variables % and ¥ and are therefore 
related to one another through the affine transformation, 
Eq. (10). According to Eq. (19), 7 controls the thick- 
ness ratio of body in axisymmetric flow, and both the 
thickness and angle of attack in plane flow. 
Moreover, the similitude based on the above analysis 
states that the properties of flow fields around bodies of 
the family defined by Eq. (19) are functions of WV ..r, 


D/P « 
/ 
M,°*” (1+ ») (G/L), 
Xx, Par Par Xia] (20) 
The corresponding equation of the shock surface is of 
the form: 


= 8[(x/L), Mar, G/L), 
P ws (21) 
Lift, moment, and drag coefficients of hypersonic 


(plane) airfoils corresponding to such bodies may also 
be expressed as functions of the foregoing parameters. 


M 
M .*Cu? = M.*k(d/L), por Por Xie] (22) 
M .*Cp 


The drag coefficient for the similar bodies of revolution 
may be written in the form 


M = f([M ar, Pos Por Xia) (23) 


where, in defining the Cp, the base area of the body is 
taken as the reference area. 
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Eqs. (19)—-(23) provide a generalization of the laws of 
similitude of Tsien and Hayes to apply for the flow 
of a real gas in thermal equilibrium about blunt-nosed 
slender or thin bodies. 

In the special case of nonreacting gas mixtures, for 
which the internal heat capacities depend only on tem- 
perature, it may be shown on the basis of Eqs. (11)- 
(16) that the restrictions for the similitude may be 
somewhat reduced, in that the specification of p ., p., 
and X;., may be replaced by that of 7., and X;,,. 
For the class of nonreacting ideal-gas flows with con- 
stant specific heats, this specification can further be 
reduced to, simply, 


Special Similitude for Blunted Wedges and Cones 


For a blunted wedge or cone, Eq. (19), governing 
the geometry of the similar bodies, is restricted to the 
form 


y/x = a (x > 0) (24) 


Here, the half cone or wedge angle a plays the role of 
7. In this case, the arbitrary length scale L in the 
transformation equation, Eq. (10), is no longer useful 
for generating new members of the family of similar 
bodies. It is convenient to choose 


L = pity) d (25) 


The functional relations constituting the similitude 
rules can now be reduced from Eqs. (20) and (21) to 
the forms 


b/P 

p/P 

S/S « M ..(y/x), M a, 

Por Per Xie) 


and 


M Pas P ws X; 0] (27) 


Note that the variables M/ ..(y/x) and M ..(6/x) can be 
replaced in the above relations, respectively, by 

It follows that the corresponding lift, moment, and 
drag coefficients of the wedges can be correlated in the 
form 


M 

M = f[M .*k(d/x), Por Por Xiwo] (28) 
M Cp) 

and the drag coefficient of the cones based on the base 
area may be correlated as 


M = f[M M Pw Po, Xia) (29) 


Approximate solutions for the problems of hyper- 
sonic flow over blunt wedges or cones have been given 
by Chernyi* for ideal gases with constant specific heats. 
The method is based on the small-perturbation theory 
and consists of introducing additional assumptions to 
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5 hypersonic flows over blunt-nosed flat plates or cylinders 
= 
- —" Examples of Correlation of Similar Hypersonic 
se Flow Fields 
Be it. n- Similitude of Flow Fields Over Flat Plates Downstream of 
05 Small Blunt Leading Edges 
Recently, Bertram and Baradell have presented in 
KH pepo . reference 31 the results of a series of very refined nu- 
= A merical analyses, based on the characteristics method, 
44: . oa 1a = for flows of an ideal gas with constant specific heats 
(y = 1.40) over a flat plate downstream of a sonic- 
— wedge leading edge.t The values of are pre- 
0 a." 10" sented in reference 31 as functions of the variable 
Mat M ..*k(d/x). This parameter was suggested by the 


Correlation of pressure increment on the surface of flat 
plate afterbodies (7 = 1.40). 


Fic. 2 


simplify the density and the pressure fields. His result 
suggests the following correlation for pressure data 
on wedges or cones: 


b'/p.M..*a* = f(é, y) (30) 


where = 


and p’ is the difference between the surface pressure 
at a ~ 0 and its corresponding value at a = 0. Ex- 
amination shows that this relation is consistent with 
the result based on Eq. (26), but is simpler. The 
similarity revealed in Eq. (30) is certainly correct in 
the flow region where the blunt nose effect is pre- 
dominant-—i.e., where £ is sufficiently large. An 
analysis based on an expansion procedure similar to 
Stewartson’s”* reveals that the error in Eq. (30) may 
be considerable in the range — < O(1). 

The laws of similitude governing the flows over the 
flat plates or cylinder afterbodies (7 ..a = 0) are also 
contained in Eqs. (26) and (27) for the blunted wedges 
or cones. Objection may arise on account of the 
ambiguity in the meaning of the affine transformation, 
Eq. (10), when tr = a = 0. This difficulty can, how- 
ever, be overcome by reducing the result desired di- 
rectly from the general similitude—1.e., from Eqs. (20) 
and (21). One can choose under this circumstance 


and a similitude will be obtained which is the same as 
that given by Eqs. (26) and (27) upon substituting 
Ma = 0. 

The solutions obtained in references 24 and 25 for 
the flat plate and cylinder afterbodies correspond to 
the case when* 


M (3+ ‘(1+y) (d x) > 1 


In view of Eqs. (26) and (27), one observes that for 
the same free-stream atmosphere (or, for the same y, 


* An error should be pointed out here in the numerical coeffi- 
cient given in reference 25 for the shape of the shock wave over 
a blunt-nosed flat plate in a flow of ideal gas with y = 1.40. The 
corrected formula reads 6/d = (0.73) k!/3(x/d)?/3, 


pressure formula of the ‘‘strong blast-wave solution” 
of references 24 and 25. As an example of the simili- 
tude discussed in the present study, this correlation is 
reproduced in Fig. 2 with the additional of a more re- 
cent set of numerical data taken from the work of 
Ferri and Pallone.** In terms of the parameter 
1/[A7..*k(d/x)], the data not only agree well with the 
strong blast-wave theory over the range where the 
theory is expected to hold; i.e., where M ..(dé/dx) ~ 
M ..[k(d/x)]'* > 1 and dé/dx < 1; but also give good 
correlation for different combinations of J/., k, d, 
etc., down to a region where 


M..(dé/dx) or (Ap/p.) < O(1) 


The theoretical basis for such correlation over the 
wide range of /.. *k(d/x) is the similitude given in Eq. 
(26). Of course, the similitude study indicates that, 
for ideal gases with constant specific heats, the pressure 
ratio on a flat plate or cylinder afterbody is determined 
only by the specific heat ratio y and the variable 
M DEY (d/x), 50 longasM ..>1,di/dx<1 
and Ap/p.. > O[(6/x)*]. The last condition results 
from the inherent error of the hypersonic small-pertur- 
bation theory previously mentioned. 

The divergence of the numerical data from the uni- 
versal trend of correlation in each series of computa- 
tions (see Fig. 2) is obviously associated with the break- 
down of the theory in the vicinity of the blunt nose. 
It is beyond the scope of the small-perturbation theory 
to account for this discrepancy. 

An attempt to correlate the coordinates of the shock 
wave downstream of the leading edge was also made 
in reference 31, and the result does not appear as satis- 
factory as for the pressure. (Readers are referred to 
Fig. 1(b) of reference 30.) It should, however, be 
pointed out that Bertram and Baradell correlated 


+ A similar computation has been given, more recently, by 


Bertram and Henderson*®® for y = 5/3 = 1.667. In the same 
reference, the calculated pressure distributions over a blunted 
plate at small angles of attack are also given for 4. = 20, 


y = 7/5 = 1.40. In view of the similitude given by Eq. (26), 
these pressure data can be used to predict or to compare with 
the pressure distributions for other values of M o, a, k, and d. 
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1.24 een can be effected even if one chooses to correlate the data’ 
fae thom in the manner of Eq. (27). 
7 eaten It should also be noted that the above inviscid blunt- 
1.20 nose effect can be realized in practice only at a suffi- 
(FT/SEC) ciently high free-stream density or Reynolds Number, 
o = ue If the density or the Reynolds Number of the free 
= Siar stream is lowered, not only the displacement effect of 
| pod the viscous boundary layer will be important in deter- 
mining the pressure field downstream of the nose, the 
“inviscid blunt-nose effect’ itself will also diminish on 
account of the interaction of the thick viscous region 
with the bow shock. Readers may refer to references 
33 and 34 for further discussions of these limitations. 
1.08 Similitude of the Real-Gas Flows Over Wedges and Cones 
1.06 : et . ct — The validity of the law of similitude extended to real 
| ® a gases can be readily checked by examining the simple 
108 ae | (BUREAU OF ORDNANCE) case of air flows over wedges. A correlation of equilib- 
= 1.107 x10 ATM. rium wedge flow data has been made in reference 22, 
i le . in which oblique shock data based on the real-gas cal- 
culation by Moeckel'‘ were correlated for the density 
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Correlations of the ‘compressibility factor’’ for 


Fic. 3(a). 
equilibrium wedge flows of argon-free air. 


5/d vs. k'?(«/d), and this correlation is appropriate 
only when the bow shock wave is extremely strong.** * 
Therefore, the correlation given in Fig. 1(b) of reference 
30 cannot be used as a test of the validity of the simili- 
tude given by Eq. (27). In fact, the shock-wave sur- 
face slopes pertaining to the shock coordinate data of 
reference 31, as close examination reveals, are generally 
so large that it is doubtful whether a good correlation 
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Fic. 3(b). Correlations of the temperature ratio for equilibrium 


wedge flows of argon-free air. 


ratio and the ratio of the shock angle to the wedge angle 
as functions of M.a at the same altitudes. In the 
present paper, similitude of the wedge flows is again 
examined on the basis of the more extensive data of 
Feldman. The ratios of temperature 72/7. and of 
density pz/p., as well as the “‘compressibility factor” 
Z» of the equilibrium air over the wedge, are accord- 
ingly correlated with respect to the parameter Ma 
under the same “‘free-stream atmosphere’’; i.e., the 
same altitude. Two typical series of the correlated 
data are presented in Fig. 3 which correspond, respec- 
tively, to altitudes of 100,000 ft. and 250,000 ft. in the 
standard atmosphere published by the Bureau of 
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Correlations of the density ratio for equilibrium wedge 
flows of argon-free air. 


Fic. 3(c). 
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Ordnance, U.S. Navy.® In order to fulfill the small- 

rturbation assumption (b), all data corresponding 
to a > 30° are excluded. Also plotted in some of these 
Figures are the correlated curves (from reference 22) for 
wedge flows of an ideal gas with y = 1.40. The small 
degree of scatter found in all correlations (Fig. 3) for 
the two altitudes examined provides a good indication 
of the validity of the laws of similitude. The correla- 
tion in pressure p/p. is not given since it can be in- 
ferred from the data of p/p., T/T, and Z already 
presented. 

One may note in this connection that the fraction of 
dissociated air over the wedge is not too substantial, as 
may be seen from the moderate deviation of the com- 
pressibility factor from unity. The corresponding 
eflects of the dissociation on the temperature and 
density, however, are by no means small because of the 
relatively large value of the heat of dissociation for 
oxygen molecules. Jn particular, the reader may note 
from these correlations that the density ratio exceeds 
the limiting value of 6 for an ideal gas with y = 1.4, 
as soon as M ..a exceeds about 6. 

Calculations of real-gas flows ove1 a cone at zero yaw 
angle have been made by Zienkiewicz,** Feldman,” 
and, recently, by Romig.** In reference 38, the ex- 
tensive data of p/p., T/T.., etc., were correlated as 
functions of M.. sin a@ for several given free-stream 
pressures P... Since it is implicit in the calculation 
that the temperature as well as the gas composition in 
the free stream remain the same throughout, these 
correlations for real-gas flows over cones are seen to 
agree with the laws of similitude given in this paper. 

The pressure over the wedge as well as the cone, 
according to references 14, 15, and 36-38, appears to be 
quite insensitive to the real-gas effect. It should be 
pointed out, nevertheless, that such a small effect in the 
pressure field is not necessarily a typical feature of 
hypersonic flows of real gases, particularly when the 
elects of surface curvature and blunt noses are taken 
into consideration. 


Concluding Remarks on the Equilibrium 
Assumption and on Simulating Nonequilibrium 
Flows 


Validity of the Equilibrium Assumption 


The assumption of local thermal equilibrium used in 
this study requires that the characteristic time for the 
activation or deactivation of the inert degrees of free- 
dom be much smaller than the flow transit time over 
the body. On the basis of studies such as references 
l7-21, a criterion for the attainment of a near-equilib- 
tum field over slender, or thin, bodies can readily be 
deduced. For example, according to the recent ex- 
perimental data of Glick and Wurster based on the 
observation of the shock thickness in oxygen (see Fig. 
2of reference 19), the relaxation time for the approach 
to the dissociative equilibrium behind the shock wave, 
after the density has been corrected to the atmospheric 
Value, appears to fall below 1 microsec. when the 


(equilibrium) temperature behind the shock becomes 
higher than 2,000°K. Accordingly, an estimate of the 
upper bound for the relaxation distance / in the flow 
field of the same gas over an airfoil, or slender body, is 


(l/U ~ 10- sec. 


where the exponent 3/2 is adopted, following refer- 
ences 17 and 21. Thus, for example, if one takes 
10, p./po ~ 10-* and ~ 10! ft./sec., the 
relaxation distance may be of the order 10 ft. or less 
and, under such circumstances, one may anticipate 
that thermal equilibrium will prevail in the flow fields 
over bodies which have a length much longer than 10 ft. 


The Possibility of Simulating Nonequilibrium Gas Flows 
Over Thin or Slender Bodies 


Finally, the possibility of simulating nonequilibrium 
gas flows over slender or thin bodies will be pointed 
out. While the arguments of Hayes® leading to the 
equivalence between the flows over a slender, or thin, 
body and a transient flow in a space of one less dimen- 
sion has been extended to include flows of real gases 
under thermal equilibrium, it is obvious that the simili- 
tude also hold among the flow fields which are “‘frozen.”’ 
In fact, inspection of the governing differential equa- 
tions for the flows of a chemically reacting two-com- 
ponent gas suggests that Hayes’ transient cross-flow 
concept may also be valid under certain conditions 
which are neither near equilibrium nor nearly frozen, 
so long as M@ > and 6 < 1. If this transient cross- 
flow concept is stipulated, a similitude for nonequilib- 
rium flows would follow, based on the possibility of 
finding a family of hypersonic flows which produces the 
same transient cross-flow field. The similitude would 
require not only that the parameters M 
(g/L), P Po, and X;., remain unchanged, but 
also that the flow transit time of the similar flow fields, 
L/U..orL/M ., be kept fixed. This appears to provide 
a way to simulate, for the same gas, and under the same 
ambient atmospheric condition, certain hypersonic in- 
viscid flows, over blunt-nosed thin or slender bodies, 
which are not in thermal equilibrium. The freedom 
of choosing the scale for the working models would, 
however, be lost. 


Appendix—Estimate of the Thickness of the 
Singular Region Near the Surface of a Blunt- 
Nosed Body and Its Effect in the Prediction of 

Pressure on the Body 


The lateral extent 6 of the singular region near the 
body surface, [see Fig. 1(c)], within which the specific 
entropy s is erroneously determined by the small- 
perturbation theory, can be estimated. The case 
chosen for examination is one of plane hypersonic 
steady flow of an ideal gas with constant specific heats 
involving very strong shock waves, and only the family 
of flows with “‘similar solutions’’ according to the small- 
perturbation theory are considered. For this case, 
the lateral extent of the shock wave envelope 6 is** 
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The lower limit of the power index m corresponds to 
the case of a flat plate afterbody. 

Let x < xy be the range of x for which the small- 
perturbation assumption fails badly is, 


for x < xy dé/dx > O(1) (A-3) 


The continuity relation gives 


from which one sees that 


bx ~ (px / px (A-5) 


where « refers to the average value (with respect to y) 
in the layer 6x. For the present purpose, one can sup- 
pose that the shock surface at x < xy can be repre- 
sented by the same equation, Eq. (A-1). It then 
follows that 6y, which is the value of 6 when dé/dx = 
O(1), can be approximately related to 6; namely, 


by/x ~ 
pw) (8% /5) ~ (8/x)™/-™ 


The thickness 5, may now be estimated by deter- 
mining the density ratio px/p., through the particle 
isentropic relation and the entropy change across the 
shock 


(A-6) 


therefore, (A-7) 


(0../px) ~ — D/(y + 1) 


The order of magnitude of the pressure on the body pz 
is the same as that behind the shock at the same station 
x;namely, 


™ Pshock/P » ~ + 1) .2(d6/dx)? (A-9) 
It follows that 


|x 
(A-8) 


pe (8/x)~° (A-10) 


Substituting Eq. (A-10) into Eq. (A-7), one has 


54/5 = (A-11) 


Because 2/3 < m < 1, the above relation indicates that 
sufficiently far downstream—-namely, for 6/x < 1— 
the layer 6 can indeed be quite thin as compared with 
the lateral extent of the cross-flow field, 6. 

As pointed out previously, the small-perturbation 
theory overestimates the shock strength and the en- 
tropy change across the shock; the incorrect values of 
entropy are inherited downstream throughout the 
layer. It is important, therefore, to see whether the 


incorrect entropy distribution near the body could 
lead to error in the solution for the pressure on the 
body. 

The y-momentum equation yields 


(A-12) 


| y=V+6% 


Ap =p 


y=Y 
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Because of Eq. (A-4), one may estimate Ap: as 


|Ap| ~ poU (bn. x) (A-13 
In view of Eq. (A-6), one has 


But the flow deflection angle near the body surface 
vp/ U’.., is always less than the shock wave angle for th, 
same Viz., 


ve/U. 0(6/x) (A-15 


Therefore, the pressure variation across the laver 4, js 


m/(l—m) 4 4 2 
Pr x 


(A-16 


which is indeed very small. Now, the pressure vari- 
ation over a distance 6x from the body according to the 
small-perturbation theory is also given by Eq. (A-16), 
For the same analysis Eqs. (A-12)—(A-16) apply. At 
the edge of the layer y = 6x(x) where the small-pertur- 
bation theory is valid, the pressure is predicted with 
an intrinsic error of order «®. It follows that the error 
in the pressure field predicted by the small-perturba- 
tion theory, for a given distance from the leading edge 
x, is uniformly of order e*, irrespective of the incorrect 
value of the entropy near the body. 
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Summary 


The problem considered is that of determining the direction 
and magnitude of the bounded thrust of a single-stage rocket, 
traveling in vacuum and a constant gravitational field, in order 
to minimize a function of the initial and final values of rocket 
mass, position and velocity components, and time. The general 
nature of the extremal arc is derived from a consideration of 
necessary conditions for the existence of a minimum “pay-off.” 
Criteria are established for the identification of extremal subarcs 
and for their composition into the complete extremal arc. It is 
found that flight takes place either at maximum or minimum 
thrust and that, at most, three such subarcs can arise. 


Symbols 


= effective exhaust velocity, constant 
= gravitational acceleration, constant 
= rocket mass 

= velocity components (Cartesian) 
= position components (Cartesian) 
= time 

= first integral, constant 

= Weierstrass function 

= G(p,4q,%, y, m, t), “pay-off” function 

= 

i 

= constant 

= —m 

= minimum 

= maximum ~ 

= real function 

= admissible variation of z 

= defined by Eq. (20) 

= Lagrangian multiplier 

= angle between thrust axis and horizontal, positive 
counterclockwise 

= constraint functions, see Eq. (3) 

orthogonal unit vectors 

unit vector along thrust axis 

vector (Ai, 

initial 

burnout 

final 

instant before corner 

instant after corner 

dot denotes differentiation with respect to time 


bin inertial frame 


(1) Introduction 


WwW: CONSIDER the problem of minimizing a func- 
tion G = Go(z0, to.) — ty), = 1,2,...0, 
subject to s independent differential constraints on the 
z, and at most 2” + 1 separated boundary conditions. 
This is the problem of Mayer as formulated by Bliss! 
with two exceptions: (1) There occur in the s differ- 
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On a Class of Variational Problems in 


Rocket Flight 


G. LEITMANN* 
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ential constraints k undifferentiated dependent vari. 
ables, called ‘‘driving functions” by Breakwell,’ where 
s<n-+k. (2) One of the driving functions is subject 
to a two-sided inequality constraint. Problems of this 
type have been discussed by Breakwell,? Garfinkel,’ 
and Miele.* 

The specific problem to be investigated is that of 
minimizing G(zZi, to, ty) where z; = p, q, x, ¥, m, 
for flight of a single-stage rocket in the vertical plane, 
Motion in vacuum and a constant gravitational field 
will be considered. There can be specified at most 2n 
+ 1 = 11 boundary conditions. Also driving functions 
8 and y cannot occur in the boundary relations. 

The equations of motion are 


p = (cB/m) cos ¥, = (o8/m)sin¥— 
and the mass flow rate, and hence the thrust, is bounded 

by 
In order to circumvent the difficulties arising from the 
fact that the thrust direction angle y is not defined 


when 6 = 0, we shall take 8, > 0, but let 6; — 0 to 
allow for unpowered (coasting) flight. 


(2) The Variational Problem 


Following the methods of Bliss' and Valentine,° Eqs. 
(1) are written 


= p — (cB/m) cos = ( 
= — (cB/m) sin y + g = 0 
=xX—p = () (3) 
= 
=. 


and Ineq. (2) is stated as 
ds = (8 — fi) (Bu — B) — y? = O (4) 


where y is a real function of time. 

There are seven physical variables, p, g, x, y, m, 8B, 
y, and five equations, [(Eqs. (3)], connecting them [in 
addition to Eq. (4), which defines y]. Thus, two de- 
grees of freedom remain, so that two variables, say 8 
and y, may be subjected to independent admissible 
variations. 

Lagrangian multipliers \;(¢), 7 = 1, 2,.. . 6, as yet 
undetermined, but assumed continuous except at a 
finite number of points of interval fo < ¢ < fy, i.e., at 
corners of the extremal arc, are introduced and the 
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fundamental function 


F Adi (5) 


i=1 


is formed. The problem is to minimize G subject to 
Eqs. (3) and (4), or equivalently to minimize 
if 
J=G+t+ f F dt (6) 
to 


where the admissible differentials, dz, and dt, at ¢ = 


and t = ty; are connected by the specified boundary 
conditions. 

The resulting Euler-Lagrange equations are 
bt+As=0; =0; 0; = 0; 
\s — (cB/m?*) (A; cos + sin YW) = 0; 
(c8/m) (Ay sin — cos = O; (7) 


(c/m) (Ay cos + Az sin — As — 
+ Bu 28) 


My = O 


Furthermore, there is the transversality condition 


dG + [Aidp + + Asdx + + Asdm — 


Cdtlii=0 (8) 
where dz, = 62; + 2,dt (9) 
C = Mp + Ag + + + Asm = constant (10) 


The first integral, C, exists, since F does not contain 
time ¢ explicitly. If function G and the boundary con- 
ditions do not contain both ¢) and ¢,, then 


(11) 


To obtain p(t), g(t), x(t), y(t), m(t), Bid), and y(t), 
corresponding to a trajectory of minimum G, Eqs. (3) 
and (7) must be integrated subject to the specified 
boundary conditions and those arising from Eq. (8). 

Some of the Euler-Lagrange equations are easily 
integrated. Thus, Eqs. (7);—4 lead to 


A3 = constant, A, = constant, 
Ai = Are — — bo) 
+ Aa(ty fj, (12) 
deo = Aa(t = te) 
Also Eq. (7)s5 yields 
tan = (13) 


Multipliers \,; and d2 are clearly continuous, so that 
tan y is continuous. However, Eq. (13) admits two 
principal solutions, 


and 


so that discontinuities in y—.e., jumps of magnitude 
™ may occur. Furthermore, there exists an ambiguity 
concerning the choice of principal solution of Eq. (13). 
These questions will be resolved by a consideration of 
further necessary conditions for the existence of a mini- 
mum value of G, namely those of Weierstrass and 
Clebsch. Before turning to these conditions let us 


(14) 
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state the relations to be satisfied at corners of the ex- 
tremal arc. 


(1) Corner Conditions 


Since p, g, 8, and ¥ do not appear in Eqs. (3) and (7), 
Dp, 9, 8, and W need be only piecewise continuous in in- 
terval ft < ¢ < ty. Thus, the extremal arc may have 
corners. At such corners the Erdmann-Weierstrass 
corner conditions must be satisfied. For the problem 
under discussion, they are 


Ae = Ke, (15) 
(16) 
Eq. (15) states that the A;,, 7 = 1, 2,... 5, have 


equal limits on both sides of a corner, and Eq. (16) 
implies the same for the first integral, C. Actually, 
the stronger condition of complete continuity can be 
seen to hold for the A;, 7 = 1, 2,... 4, in view of Eq. 
(12). 


(4) Weierstrass Condition 
A necessary condition? for the existence of a mini- 
mum value of G is 
(17) 


where 
E F(z,*, 3,*) F(2;, Zi) 

8 

— (OF/28,) (18) 
and * denotes functions subjected to finite, admissible 
variations. The z,;, the derivatives of which occur in 
F (here p, g, x, y, m), must be continuous—i.e., p* = 
= q,x* = x, y* = y, m* = m. 

Evaluation of Ineq. (17) leads to 


Bx — B*x* > 0 (19) 
where 
k = (c/m) (A, + sin — As (20) 
x* = (c/m) (A, cos y* + sin — As! 


Ineq. (19) must be satisfied for all admissible vari- 
Thus, for 


pry 
Ineq. (19) becomes 
cos + sin > A, cos + de sin (22) 


ations, not all zero. 
(21) 


Likewise, for 


B~p*, p=y* (23) 
Ineq. (19) is (8 — B*)x«>0 (24) 
For 8 = fu, (A: — B*)« >0 (25) 
and for 8 = Bu, (8. — B*)« > 0 (26) 


¢ See also condition (11) of reference 2. 
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We can now draw a number of important conclu- w= 0,7 #0, ie, B<B< By; 
sions. In order to satisfy Ineq. (22) for all admissible Ao ¥ 0,7 = 0, ie, B= 6: or B= Bu; (38) 
variations in must be chosen so as to maximize A= 0,7=0, te, or 
Ai cos YW + sin Hence, the extremal arc consists of subarcs of minimum, whe 
with respect to y. It is seen that this value of y is ond apes has 
f woe been mentioned by Miele* for flight in a_ resisting 
precisely one of the principal values of Eq. (13). Fur- : . It 
thermore, the choice of principal value of Eq. (13) is : : and 7 
ig ; We shall now show that either (a) the solution js 
now unambiguous—namely, it is the one which makes Sine 
nonunique, including the case for which the value of Gf 4 
Ai cos YW + Ax sin y > 0 (27) does not depend on 8(¢); or (b) the solution is unique but | °°” 
case: 
Since B=, or B= By tion 
B*< (28) and never 6; < < B,; i.e., only subares of minimum or | 
axi thrust can arise. This corresponds to the 
Ineqs. (25) and (26) show that 
(35) findings of Lawden® and Okhotsimskii and Eneev,’ for 
k <0 = problems with unbounded thrust. We shall show, 
and x 20 when 8 = 8, sn too, that for Case (b) the extremal arc can be com- 
posed of no more than three subares. 
(5) Clebsch Condition Upon differentiating « and invoking Eq. (7),-5 it is} A"? 
seen that sugs 
Another necessary condition for the existence of a i= + (39) 
minimum value of G is 
8 8 If = A» = O and « = ina nonzero interval, then and 
(oF > (30) A, = 0,7 = 1, 2,...6, everywhere. 
iti If not both A; = 0 and Ay = O, then if 
where the z, = p, g, x, ¥, m, 8, y, y, and the z,, not all a (40) whet 
zero, constitute a set of values satisfying the con- | Ther 
straints Eq. (39) yields tan Y = —(A3/Aq) (41) 
8 E D2): ‘ the 
‘ : : q. (41) together with Eqs. (12) and (13) show that : 
(084/081), (31) there exists at most one instant of time whenx’ = 0). 
tion. 
In view of Eq. (7), Ineq. (30) is poin' 
k constant (42) and. 
(cB/m) (Ay cos + d2 sin ¥) — exce} 
+ (,)?] > 0 (32) and, in particular, 
For m= 7, = 0, and #0 (33) k #0 (43) of 
tinui 
Ineq. (32) states that Thus, in view of Eq. (7)z, Eqs. 
A, cos + Ao sin y > O (34) As 0) (44) 
which is the same as Ineq. (27). For and Bi + Bu ¥ 28 (45) | so th 
tm or wr, #0, and mr = 0 (35) Consequently, only condition (38). can arise; only } at sc 
: subares of minimum or maximum thrust can occur. 
Ineq. (32) yields If 
so th 
Ae < (36) As = = O (46) 
However, from Eq. (7)z, Ai; = constant, A» = constant (47) 
A si 
= — 9 a7 
K Ae(Bi1 > Bu 26) (37) so that y = constant (48) App 
Upon putting 6 = 6, and 8 = £, in Eq. (37), it is seen Thus a (49) that 
that Ineq. (36) again leads to conditions (29). ‘ ag , As 
and kK = constant (50) that 
(b) k 
(6) Composition of the Extremal Arc If « ¥ 0, condition (38), is again the only possible one. two 
In order to satisfy Eq. (7)s one of the following condi- If TI 
tions must hold: (51) Ore 
+ It may be noted that these conclusions hold also for 8; = all of conditions (38) could arise. However, now Eqs. "7 
0, as do the results of Section (5). (7)5 and (7); become, respectively, retait 
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+ (m/m)K = 0 (52) 
Asm = K (53) 

where 
K = c(A\, cos + Sin = constant (54) 


It is seen that Eqs. (7)5 and (7)7 are not independent 


~ and that an arbitrary m(t) satisfies this set of equations. 


Since m(t) is arbitrary, A(t) is arbitrary, save for the 
condition of being bounded. Consequently, for those 
cases when \; and )\y both vanish and x = 0), the solu- 
tion is nonunique,} including cases for which G does 
not depend on @(t). Some simple and hence trans- 
parent examples of such behavior are discussed in the 


Appendix. 
As stated above, y must be chosen so as to maximize 


cos + Ae sin 


An interesting interpretation of this condition has been 
suggested by J. V. Breakwell. Consider the vector 


A= Ai + Ay (55) 
and the unit vector along the thrust axis 
u=icosy+Jjsiny (56) 


where i and j are a pair of orthogonal unit vectors. 
Then 


= A, cos + Sin (57) 


which is maximum when ) and u point in the same direc- 
tion. Hence, w is the direction angle of \—i.e., A 
points along the instantaneous thrust axis. Since ), 
and \» are continuous, y is continuous with one possible 
exception—namely, if and pass through zero simul- 
taneously. At such an instant y experiences a jump 
of magnitude 7, x is continuous, but « has a discon- 
tinuity consisting of a sign change. We note that 
Eqs. (12) yield 


— Avo = (Ar — Ato) (58) 


so that if A = de = () (59) 


at some instant, then 
doo/Aro = Aa/As (60) 


so that 
tan Y = A4/A3 = constant (61) 


A simple example of such a case is included in the 
Appendix. In view of Eqs. (39) and (61) it is clear 
that x ¥ 0 for the case \; = Ay = 0 at some instant. 

As a result of the preceding arguments it is concluded 
that if (a) « is continuous, « can have at most two zeros; 
(b) x has a discontinuity, x # 0, and « can have at most 
two zeros. 

Thus, in either case, the extremal are consists of no 
more than three subarcs: 


+ This difficulty does not arise if first-order terms in g(x, y) are 
retained. 


By Bu = Bi or Bu By Bu, or 
8 or B&—f& or or 


A consideration of the corner conditions leads to a fur- 
ther result. Eqs. (16) and (7)s yield 


= Bxi+ (62) 
where =f, B+ = B, 
or B+ = 
However, K—- = K+ (63) 


whence it follows that 
= = (64) 


at corners between subares. We note that burn-out, 
m = my, must occur either (a) at a corner preceded by 
a powered arc (8 = @,) and followed by the final coast- 
ing are (8 = 6,0), and ¢t; > ¢,; or (b) during the final 
powered are (8 = @,), so that ¢; = % This can also 
be seen from the following argument. If /¢,;> ¢,, a final 
coasting are is required. However, if burn-out, ¢ = 
t,, occurs before a corner—-i.e., when « > 0-—then no 
coasting are can follow, since « < ( during coasting. 
Hence, burn-out must take place at x = () or during a 
final powered arc. 

Furthermore, if 4; = A» = 0 at some instant, then 
k = constant when x < 0). 

The following possibilities (Figs. 1 through 6) can 
arise. 


B 
t . t 
Fic. 1 Case A. 
Case A: 
(a) my = my, (coasting flight only), tf; < 4; 
(b) mo > my, ty = ty > th. 
« B 
H 
+ t 
Fic. 2. Case B. 
Case B: 


(a) m > my, ty = ty Sh: 
(b) my > my, te = 


(41) 

y that = 
= 0), 
(42) 

(44) 
ur. 

(47) 
(48) 

(50) 

one. 

(51) 
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Fic. 3. Case C. 


Case C: 


(a) my = my (coasting flight only), ty < 4; 


(b) My Po My, ty > to, th = to; 


(c) my > My, h< ty < bo, t, = ty. 


« B 
WA 
+ 
+ to 
Fic. 4. Case D. 
Case D: 
(a) my > My, ty = ty hh; 
(b) My > Ny, ty > te, th = ty; 
(c) mo > myth < ty < te, ty = hh. 
B 
+ t, 
Fic. 5. Case E. 
Case E: 
(a) me = my (coasting 
flight only), ty < t, 
(b) mo > my, 
= »o = ( = 
ty > ty ts de ) at i 
(c) mo > my, 
ts, = ty 
B — 
' 
+ 


Fic. 6. Case F. 


1959 


Case F: 


(a) m >m,tr=t 

(b) mo > my, ty > ts ty = ty 

(c) mo > my, 
=h 


(7) Conclusions 


The following principal conclusions have been de. 
rived concerning optimum trajectories of single-stage 
rockets in vacuum, constant gravity flight: 

(1) The optimum trajectory is composed of maxi- 
mum and minimum (coasting) thrust arcs only.t 

(2) There occur at most three subarcs in a trajectory, 

(3) An arc is one of maximum thrust when « > 0, and 
of minimum thrust (coasting) when « < 0. 

(4) Transition between arcs takes place when x = 

(5) Burnout (m = my) occurs either at the junction 
between a thrust and final coasting arc, or during a 
final thrust arc. 


Appendix—Examples 


Here we shall discuss briefly some very simple 
examples merely with a view toward making some of 
our conclusions ‘‘palatable.” 


Example (1): Maximum Horizontal Velocity 


G = 


= po 
q = d, t=t,{m =m, 


m = Mo 
In this case, Ae = Ae ='0 
so that Ai = constant = Aj; = 1 
A2 = constant = Ao, = O 


Consequently, as expected, 


tany = 0, and y=0 
Also, <= 0, and C= px 
If ¢; is unspecified, C=0 
so that k = 0 


According to the results of Section (6) this implies that 


G = —p; does not depend on @(t). Since 
ps = cln + po, = 0) 
7 Nore: After submission of this paper a report by Miele 


and Capellari® came to the author’s attention. Of interest are the 
examples of the problems of attaining maximum range and maxi- 
mum altitude. For these problems the nonexistence of variable 
thrust arcs is shown and the probability of their nonexistence in 
other cases surmised—a surmise which was indeed correct as 
shown here. 

t An equivalent statement for a particular problem has been 
found previously—e.g., see Miele.’ 
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this is clearly the case. As a matter of fact, the result 
holds whether or not t; > t, is given (i.e., C = 0). 


Example (2): Maximum Vertical Velocity 


Pp = po 
q = % 
t=t4=0 x =X, t= t,{m = my, 
y = Vo 
m= mM 
In this case, hes SO 
so that Ai = constant = Ay = OV 
= constant = As, = 1 


Hence, as expected, 


tany = ~, and y= (1/2)r 
Here again, t= 0 
Also, C= 
lf ty is unspecified, 
C=0, and =g 
Le., k = constant > 0 


Hence, the trajectory consists of a maximum thrust arc 
only. In other words, ty = ¢,, where f is the shortest 
burning time compatible with 8 < @,. Again this is 
obvious from 


= (1/2)m] 


For the case of unbounded thrust, tr = 4 = 0 (im- 
pulsive boost), whereas here 


= €ln (mo/my) — gty + qo, 


ty = (mo — my)/Bu 


this fact should be predicted by the results of Section 
(6). It can be shown, e.g., see Breakwell,” that 


If ty; > ty is specified, is independent of @(¢) and 


C = 
but Oq,/Oty = —g 
so that Bx = 0 
and, hence, xk = 0 


This, however, is precisely the criterion for G not de- 
pendent on £(t). 


Example (3): Maximum Payload 


G my 
= po 
q = 
m= mM 
In this case, 
= Axe = 0, Asy = 1 


tan = = constant 


This is a case for which x possesses a cusp, here at ¢ = 
ty, and, hence, either Fig. 5 or Fig. 6 of Section (6) 
applies. Since 


kp = = -1 <0 


the final arc is a coasting one, so that Fig. 6 portrays 
the situation. 
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Summary 


A general method is presented for investigating optimum con- 
ditions of the quasi-steady mechanics of flight. 

For motion in a vertical plane the problem of extremizing an 
arbitrarily specified function of altitude, Mach Number, lift, 
path inclination, engine control parameter, and thrust inclination 
is considered for an aircraft which has to satisfy the equations of 
motion and 3 additional arbitrary constraints. 

A generalized solution is obtained in a determinantal form. 
The characteristic of this solution is that it unifies into a single 
equation the results of a large segment of the previous contri- 
butions to the quasi-steady mechanics of flight. Particular 
problems such as maximum speed, maximum range, maximum 
endurance, ceiling, steepest ascent, best rate of climb, flattest 
descent, ete., are thus all covered by the same determinantal 
equation. 

The optimum ratio (RX) of induced drag to zero-lift drag is 
evaluated for arbitrary relationships between zero-lift drag 
coefficient, induced drag coefficient, thrust, specific fuel consump- 
tion, and Mach Number. 

Design problems are also investigated, such as those associated 
with the selection of the best wing surface or the best aspect 
ratio for a given performance to be optimized. 

Finally, the paper is completed with the study of the optimum 
flight conditions for curvilinear motion in a horizontal plane. 
Another general determinantal equation is derived, from which, 
as a particular case, the turning flight with maximum angle of 
bank, maximum angular velocity, or minimum radius of curvature 
is investigated. 


(1) Introduction 
A’ IMPORTANT BRANCH of the aeronautical sciences 
is the mechanics of flight of an aircraft. Its 
overall scope is to determine the performance of a fly- 
ing machine; its particular objective is to detect those 
special conditions under which an optimum is reached 
with regard to both design and flight operations. 
It was a common belief among engineers, in the years 
preceding World War II, that the mechanics of flight 
had reached a conclusive stage of development. This 
opinion has been proved incorrect. In fact, the ad- 
vent of jet engines as aircraft propulsion systems and 
the increase in flight velocities have generated a new 
wealth of important and unsolved problems of applied 
mathematics. This concept particularly applies to 
rocket-powered missiles. These missiles, because of 
the eminently non-steady character of their motion, 
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have originated an ever-increasing demand for a 
variational approach to the mechanics of flight. 

Nevertheless, with regard to aircraft propelled by 
air-breathing engines in general (and, in particular, 
with regard to turbojet powered aircraft, which will 
dominate the next decade of commercial air transporta- 
tion) there are still important flight conditions for which 
the acceleration terms may be neglected. These flight 
conditions are consequently amenable to a less sophisti- 
cated, but useful, treatment by the ordinary theory of 
maxima and minima.{ 

In the present paper, the quasi-steady point of view 
is adopted and a re-analysis is offered of the problems 
of the mechanics of flight at large. Several simul- 
taneous objectives are embodied in this investigation: 

(a) To formulate the optimum conditions of flight 
in a general fashion, valid (among other factors) for an 
arbitrary dependence of the characteristics of the air- 
craft and of the engine on the Mach Number. 

(b) To unify all the optimum conditions of the me- 
chanics of flight in a vertical plane into one single 
equation; this equation contains, as a particular case, 
the answer to a variety of problems such as ceiling, best 
range, best endurance, maximum level speed, steepest 
ascent, best rate of climb, flattest descent, etc. 

(c) To unify all the optimum conditions of the me- 
chanics of flight in a horizontal plane into another 
single equation, analogous to the one discussed in (b). 

(d) To approach design problems such as_ those 
associated with the selection of the best wing surface 
or the best aspect ratio for a given performance to be 
optimized. 


(2) Fundamental Equations 


Denote with 7 the thrust, D the drag, L the lift, W 
the weight, 6 the inclination of the velocity with respect 
to a horizontal plane, and w the inclination of the thrust 
with respect to the velocity. Denote with t = p/ps 
the ratio of static pressure at the altitude / to static 
pressure at the altitude 4, of the tropopause. Indicate 
with V the flight velocity, a() the speed of sound, and 
M = V/athe Mach Number. Assume a drag function 
of the form D = D(z, M, L) and a thrust function T = 


t Incidentally, this writer has recently shown? that— 
when the acceleration terms are neglected—the results obtained 
via calculus of variations are identical to those which can be 
obtained with the ordinary theory of maxima and minima. 
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\/, 8), where is a variable controlling the enginef 
performance and having a function analogous to the one 
of the accelerator pedal in an automobile. 

With the previous considerations in mind, the funda- 
mental equations of the quasi-steady flight along the 
tangent and the normal to the flight path are written as 


follows: 


VM, L, 6, 8, w) = M, B)cos w — 


D(x, M, L) Wsing@=0 (1) 
L, 0, 8B, w) = M, sin w + 
L— Wecos@=0 (2) 


For a given weight W, the above two algebraic equa- 
tions involve the six variables 7, M, L, 0, 8, and w. 
Four degrees of freedom are left and, as a consequence, 
some optimum requirement can be imposed on the set 
of solutions of Eqs. (1) and (2). 


(2.1) Additional Constraints 


In many engineering applications it is of paramount 
interest to study particular solutions of the algebraic 
system formed by Eqs. (1) and (2), more specifically, 
those solutions which simultaneously satisfy a set of 
additional constraints.~ The latter are symbolically 
indicated asfT 


®;(7, M, 6, B, w) = 0) 
P,(7, M, L, 6, B, w) 0 
;(7, M, L, 6, B, w) 0 


(3) 


After denoting with m the number of degrees of free- 
dom associated with Eqs. (1) to (3), the following 
situation arises: m = 1, if none of the additional con- 
straining functions is identically zero; n = 2, if only 
$; is identically zero; n = 3, if both # and ®; are 
identically zero; and n = 4, if ®3, &4, and ®; are identt- 
cally zero. 


(3) Statement of the Problem 


A functional expression of the form 


Vv = V(n, M, L, 0, B, w) (4) 


is now considered and the minimal problem formulated 


as follows: ‘‘Among all sets of solutions of Eqs. (1) 


t Within the context of the present paper the variable 8 is 
termed engine control parameter, thrust control parameter, or power 
setting. Nevertheless, its actual physical meaning is not speci- 
fied, because of the great variety of existing types of air-breathing 
engines and in order to develop a general theory. Asan example, 
however, 8 can be identified with the number of revolutions 
of the turbine-compressor group for a turbojet engine having 
fixed geometry or with the position of the fuel control lever, etc. 

t Within the context of the present paper the constraints (1) 
and (2) are termed fundamental; on the other hand, the con- 
straints (3) are termed edditional. 

tt As an example, a level flight condition is expressed by 6 = 0. 
As another example, a condition in which the thrust control 
parameter is specified is represented by 6 — const. = 0. 


to (3), find that particular set which extremizes (i.e., 
maximizes or minimizes) the V-function.’’ 

To this effect, the constant Lagrange multipliers 
Ai. . . Ag are introduced and the following composite 
expression formed: 


(5) 


where Ay = 1 and &;=VW. It is known from the theory 
of the Lagrange multipliers that the necessary conditions 
for an extremum are 

OF/dz; = 0 (j= 1,..., 6) (6) 
where 2; = 7, 2. = M,23 = L, % = 0,2 = 8, and ss = w. 
The development of Eq. (6) leads to 


[(O7/Om) cos w — (OD/Om)| + A(OT/Om) sin w + 


(0/07) (A3®3 AgPy + As®s (7) 

[(O7/0M) cos w — (OD/OM)| + 
sin + 

(0/OM) (As®3 + + + V) = O (8) 
—Ai(OD/OL) + Ax + 

(O/OL) (As®3 + + + = O (9) 
— ),W cos 6 + sin + 

(0/086) (A3P3 + + As®s = 0 (10) 
\(O7T/0B8) cos w + A(OT/OB) sin w + 

(0/OB) (AsPs + + ASD, + V) = O (11) 
sin w + cos w + 

(0/Ow) (AzP3 + AsPy + ASPs + Vv) (12) 


(3.1) Explicit Form for the Optimizing Condition 


In order to formulate the optimizing condition in an 
explicit fashion the 6 X 6 matrix 


is considered [i = 1,...,6;j7 = 1,..., 6). Its ele- 
ments A, are the partial derivatives (with respect to 
the variables of the problem) of the first members of 
the constraining equations &,, ®;) and of the 
function to be extremized (®,= W). 

(3.1.1) Problems With One Degree of Freedom— 
Assume now that the first members of the three addi- 
tional constraints are not identically zero (nm = 1). 
After eliminating the Lagrange multipliers from Eqs. 
(7) to (12) the optimizing condition is explicitly stated 
as follows: ‘“‘the determinant of order 6 associated with 
the square matrix (13) 1s to be zero’’; 7.e., 


(13) 


tt In practical cases the ¥-function can be the flight velocity, 
the flight altitude, the path inclination, the rate of climb, the 
range per unit fuel consumed, the endurance per unit fuel con- 
sumed, etc. 
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This equation, in combination with the fundamental 
constraints (1) and (2) and with the additional con- 
straints (3), completely determine the solution of the 
optimum problem. 

(3.1.2) Problems With Several Degrees of Freedom— 
If the problem admits degrees of freedom, the elimi- 
nation of the Lagrange multipliers from Eqs. (7) to (12) 
leads to the following general rule for the optimizing 
condition: ‘“‘extract from matrix (13) any n non-iden- 
tically zero determinants of order 7 — n; set the above n 
determinants equal to zero.”’ 


(4) Analysis of Some Special Problems 


The dominant characteristic of the previous solu- 
tions is that they unify all the optimum conditions of 
the mechanics of flight in a vertical plane into a single 
generalized treatment. The answer to_ particular 
problems is simply found by properly specializing the 
four functions ®;, ®,, ®;, and VY. This concept is illus- 
trated in the following sections. 


(4.1) Optimum Inclination of the Thrust Axis 


It is desired to maximize the speed in level flight for 
a given altitude and a given thrust control parameter. 
After setting 


= M a(x) 
= 0 
®, = 7 — const. = 0 (15) 
&, = B — const. = 0 
the determinantal equation (14) becomes 
ie oF sin w 0 1 0 
us Tv Tv 
—sinw 000 a 
a 1 000 0 =0 
OL 
0 100 0 
oF cos w —sinw O00 1 0 
—T sin w Tcosw 000 O (15a) 


oT OD oT . Of; OF, OF; OV 

-cosw—- — —snw — — — — 

OT Or On Or OF 

OT Oty Ob av 

OM OM OM OM OM OM OM 

OL OL OL OL OL 

—W cos 6 W sin @ > = 

08 O08 08 

Of; Of; Of; OW 

—— cos w sino — — — — 

OB OB OB 

—T sin w T cos w 
Ow Ow Ow Ow | (14) 


1959 


admitting, therefore, the following solution: 
w = arctan (OD/OdL) (16) 
A drag polar of the form 
Cp = Cpn,(M) + K(M)C;? (17) 


is now considered, where Cp is the total drag coefficient, 
C, the lift coefficient, Cp, the zero-lift drag coefficient, 
and K the induced drag factor. After accounting for 
the classical definitions of lift and drag, Eqs. (16) and 
(17) lead to 


w = arctan (2KC;_) (18) 
in general and to 
w & 2KC, (19) 


for small lift coefficients. The following remarks are 
in order. 

(a) The solution expressed by Eq. (19) was origi- 
nally derived by Riparbelli! for the particular case of a 
piston-engined aircraft with Cp, and A independent 
of the Mach Number. Actually, Eq. (19) holds, what- 
ever be the power plant and the relationship between 
the zero-lift drag coefficient, the induced drag coeffi- 
cient, and the Mach Number. 

(b) The optimum inclination of the thrust axis [Eq. 
(16)] holds not only for level flight, but in general for 
all those flight conditions where the additional con- 
straints and the function to be extremized are such that 


0®,/OL = 0, 04,/dw = 0 (i = 3, 4, 5, 6) (19a) 


(c) The optimum inclination of the thrust axis in 
unaccelerated flight [Eq. (16)] is formally identical 
with the analytical expression which can be calculated 
by more sophisticated variational procedures for 
accelerated flight.” * 


(4.2) Maximum Range at a Given Altitude 


Consider an aircraft operating in level flight at a cer- 
tain prescribed altitude} and assume that the thrust 
and the velocity are parallel. Denote with c(z, M, 8) 
the specific fuel consumption—i.e., the weight of fuel 
consumed per unit time and unit thrust. 

The problem of the maximum instantaneous range 
per unit fuel consumed consists of extremizing the 
quantity 


W = Ma(nm)/c(x, M, B)T(a, M, B) (20) 


subject to the equations of motion and to the following 
additional constraints 


®; = 7 — const. = 0 
= 0 (21) 
= w 


One degree of freedom is left and the optimizing con- 
dition is consequently expressed by Eq. (14), that is, by 


T The engine control parameter 8 is not prescribed in the present 
problem. 
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The above expression leads to 
| &T — D)/oM (0/0M)(Ma/cT) | 
| = 0 (22) 
oT /op (0/08)(Ma/cT) | 


arelation which, after accounting for the constraining 
equations, yields 


(0 log c/08) [0 log (T/D)/0 log M] + 
(0 log 7/08) {1 — [d log (cD)/d log M]} = 0 (23) 
The ideal case of an engine such that 0c/d06 = 0, 


07/08 ¥ 0 is now considered. Eq. (23) consequently 


reduces to 
0 log (cD)/0 log M = 1 (24) 


For a polar obeying Eq. (17) the logarithmic derivative 
of the drag with respect to Mach Number takes the 


form 

dlogD 

Ylog 

2+ (d log Cp,/d log M) + R[-—2 + (d log K/d log M)] 


1+R 
(25) 
where R = K C,?/Cp, is the ratio of induced drag co- 


dlicient to zero-lift drag coefficient. As a consequence, 
Bq. (24) is transformed into 


log (cCp,)/0 log M] 


= 2 
3 — [0 log (cK)/0 log M] 


Notice that the optimum Mach Number is to be con- 
sistent with the relationship 


/ 
We 1 + [0 log (cCp,)/0 log M] (27) 


K 3 — [0 log (cK)/0 log M] 


where y is the ratio of specific heat at constant pressure 
to specific heat at constant volume (air), p the local 
atmospheric pressure, and S a reference surface for the 
aerodynamic coefficients. Consider now a long-range 
ttuise trajectory and assume that the optimizing condi- 
tion (27) is maintained at all time instants. The follow- 
ing situation arises: owing to the consumption of fuel, 
the weight W varies along the path. As a consequence, 
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the optimum Mach Number varies, generally decreas- 
ing as the flight progresses. 

When the Mach Number derivatives of c, K, and 
Cp, are ideally zero, Eqs. (26) and (27) reduce to 


R=1/3 


(28) 
(29) 


i.e., to the conditions indicated, among others, by Ash- 


. kenas,* Santangelo,’ and Lippisch® in connection with 


the flight performances at relatively low speed. 

Remark—For a piston-engined aircraft the specific 
fuel consumption referred to the unit thrust (c) and the 
specific fuel consumption referred to the unit power 
(cp) satisfy the relationship 


c= Cp(Ma/n) 
where 7 is the efficiency of the propeller. Assuming that 


the ratio cp/n is ideally independent of the Mach Num- 
ber, one obtains 


0 log c/0 log M = 1 
Consequently, the optimum ratio of induced drag to 
zero-lift drag becomes 
R=1 


in the region where the derivatives of K and Cp, with 
respect to Mach Number are zero. Furthermore, the 
optimum Number becomes 


M = WK/Cp, V2W/ypS 


(4.3) Maximum Endurance at a Given Altitude 


The problem of the maximum instantaneous endur- 
ance per unit fuel consumed consists of extremizing 


WV = 1/c(z, M, B)T(z, M, B) (30) 


subject to the constraints (1), (2), and (21). The num- 
ber of degrees of freedom being = 1, Eq. (14) reduces 
to 


o(T — D)/oM 
oT 


If the specific fuel consumption in independent of the 
thrust control parameter, then Eq. (31) yields 


0 log (cD)/0 log W = 0 


For a parabolic drag polar, the following results are 
obtained: 


(0/OM)(1/cT) | 


=0 (31) 
(0/08) (1/cT) 


(32) 


2 + [0 log (cCp,)/O0 log 


33 
2 — [0 log (cK)/0 log M} (33) 


tT If the overall range is desired, the following differential equa- 
tion is to be integrated 


where the instantaneous value of WV results from the solution of 
the Jocal optimum problem. The integration process leads to 
closed form expressions only for particular cases, for example, 
the case represented by Eqs. (28) and (29). 


| 
M = W3(K/Cp,) V2W/ypS 
= 
a cer- 
hrust 
M, 8) 
range 
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7 
(21) 
con- 
s, by dw 
’ 
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(34) 


20 
ypS [0 log (cK) log M| 


When the Mach Number derivatives of c, Cp,, and K 
are ideally zero, the previous expressions lead to (see, for 
instance, reference 5) 


R=1, M=WK/Cp,V2W/ypS (35) 


(4.4) Maximum Range for a Given Thrust Control 
Parameter 

Consider the problem of maximizing the expression 
(20) for the case where the engine control parameter is 
prescribed.t Assuming three additional constraints 
of the form 


= — const. = 


= 0 (36) 
=w = 0 
the op imizing condition [Eq. (14)] reduces to 
| — D)/dr (0/dm)(Ma/cT) 
| — D)/OM (0/0M)(Ma/cT) 


(4.4.1) Turbojet Aircraft Operating in the Strato- 
sphere—Attention is now focused on an ideal turbojet 
engine operating in the isothermal region of the strato- 
sphere (a = constant) according to the laws 


M), Cx (8, M) 


where the subscript « denotes magnitude evaluated at 
the tropopause. 

For a parabolic polar, the logarithmic derivative of 
the drag with respect to the pressure ratio takes the 
form 


(38) 


0 los D/O log = (1 — R)/(1 + R) (39) 


As a consequence, Eq. (37) yields the optimum ratiof 
of induced drag to zero-lift drag’ 
_ 2+ fo log (Cp,/T )/O log M | 


4 — [0 log log M] 


(40) 


In turn, the Mach Number is to be consistent with the 
fundamental equationff (see reference 7) 


2T x 6 — [0 log (Kcx?T%?/Cp,)/0 log M] 


—— =Cp,M? 
Co, 4 — [0 log (Kcx?T%)/0 log 


(41) 


7 The flight altitude (i.e., the pressure ratio +) is not prescribed 
in the present problem. 

¢ The optimum ratio of induced drag to zero-lift drag can be 
considerably different from the value obtained from the low speed 
flight theory. In a typical case, R may range from 1/4 to 4 
times the subsonic value. 

tt For a given thrust control parameter, Eq. (41) may admit 
either one or several solutions—i.e., several Mach Numbers 
yielding a stationary range—depending upon the thrust level 
and the aerodynamic characteristics of the aircraft. When 
several solutions exist, a comparison of their relative merits 
is to be carried out by evaluating the W-function, supplied by 
Eq. (20). 


THE AERO/SPACE 


SCIENCES—SEPTEMBER, 


1959 


Consider now a long-range cruise trajectory and as. 
sume that the optimizing condition (41) is maintained 
at all points. In view of the fact that the instantaneous 
weight W does not appear in Eq. (41), one deduces that 
the cruise is performed at constant Mach Number (/f 
and at constant ratio (R) of induced drag to zero-lift 
drag. The instantaneous lift coefficient C, being g 
constant, it follows that the ratio W/7z has a constant 
value as the flight progresses; this concept is the origin 
of the cruise-climb technique already outlined by 
Ashkenas‘ for the case where Mach Number effects on 
K, Co,, Cx, Tx are not considered. Rigorously speaking, 
the above cruise-climb contradicts the hypothesis 6 = (), 
Actually, however, the values of @ are so small (order 
of magnitude: 10~* radians) that the weight com- 
ponent on the tangent to the flight path is less than | 
per cent of the thrust. The assumption 6 = 0 is, there- 
fore, more than justified with regard to practical engi- 
neering applications. 

For the ideal case where the Mach Number deriy- 
atives of Cp,, K, cx, and 7, are zero, Eqs. (40) and 
(41) lead to® ® 


R= 1/2, M = (42) 


(4.5) Maximum Endurance for a Given Thrust Control 


Parameter 


For the problem of maximizing the quantity (30), 
subject to the equations of motion and to the additional 
constraints (36), the optimizing condition is given by 

| — D)/dm (0/Om)(1/cT) 
=() (43) 
o(T — D)/oM (0/0M)(1/cT) 
For a parabolic drag polar and for a turbojet engine 


operating in the isothermal stratosphere according to 
Eqs. (38), the optimizing condition leads to 


2 + [0 log (Cp,/Tx)/d log M} 


(44) 
2 — [0 log (Kcex?7T%)/0 log M| 
2 — [0 log log M] 


(45) 

When the Mach Number derivatives of Cp,, K, 7%, & 
are zero, the above expressions reduce to 

R=1, M = VTx/vpxSCo, (46) 


implying that the optimum operating altitude is 
identical with the theoretical ceiling of the aircraft.{t 


(4.6) Altitude for Maximum Speed in Level Flight 


Assume now that the thrust and the velocity are 
parallel, that the trajectory is horizontal, and that the 
engine control parameter is given. The altitude for 
which the speed is stationary is to be found by extrem- 
izing 


tt The application of the present general method to the study 
of the ceiling or the fastest ascent” is omitted from the present 
article for the sake of brevity. 
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WV = (47) 


subject to the constraints (36) and to the equations of 
motion. The development of the determinant (14) 
leads to the important equation 


(0 log (7D) log x] — 
(d log a/d log m) [0 log (7, D)/0 log M| = 0 (48) 


This equation has no solution for a turbojet aircraft 
operating in the tsothermal stratospheret in accordance 
with the hypotheses of Section (4.4.1). Consequently, 
two alternative possibilities must be inferred—either 
the altitude for maximum speed is located in the tropo- 
sphere or the highest value for the velocity (not a maxi- 
mum in the analytical sense) is reached at the tropo- 
pause. 


(4.7) Steepest Climb 


The unaccelerated steepest climb problem consists 
of extremizing 


(49) 


subject to the constraints 


®; = 6 — const. 0) 
?, = 7 — const. = 0 (50) 
=w = 0 


and to the equations of motion. The determinant 


(14) is once more applied, leading to 
0 log (7°, D)/0 log M = 0 (51) 


It should be noticed that this result also holds for the 
flattest descent of a glider (7 = 0), leading to 


0 log D/O log WJ = 0 (52) 
in general and to 


2+ (d log (53) 
2 — (d log K/d log ) 

for a parabolic drag polar. Consider now the all- 
supersonic domain and assume, for simplicity, that 
Cp, is proportional to /*-° (where x is a constant) 
while K is proportional to 1/7. Eq. (53) yields R = x, 
implying that the induced drag for flattest glide is equal 
to the zero-lift drag for the linear drag law (x = 1), but 
twice the zero-lift drag for the quadratic drag law 


fe = 2). 


(5) Design Considerations 


The method previously outlined yields the best 
flight performances for a given aircraft configuration 
and for a given engine. The same method can be em- 
ployed to solve important design questions. 

As an example, consider the problem of determining 


+t Numerical analyses show that the speed for level flight is a 
monotonically decreasing function of the altitude. 


the optimum wing area (.S), i.e., the wing surface which 
extremizes the function (4), subject to the constraints 
(3) and the equations of motion.{ Two new elements 
must be considered—namely, a drag function D = 
D (m, M, L, S) and a weight function W = W(S). As 
a consequence, the set of six Eqs. (7) to (12) must be 
completed by a seventh equation, relative to the opti- 
mization of S. After assuming ®,; (7 = 3, 4, 5, 6) inde- 
pendent of .S, this seventh equation reads as follows: 


\,[(OD/OS) + (dW/dS) sin 6] + 
he(dW/dS) cos6 = 0 (54) 


For 0%,/0L = 0 (i = 3, 4, 5, 6), Eq. (9) simplifies into 
\(OD/OL) = de = 0 (55) 


Elimination of the multipliers from Eqs. (54) and (55) 
yields the important relation 


(OD/OS) + (dW/dS) [sin @ + (OD/OL) cos 6] = 0 
(56) 


Despite the generality of the above result, a word of 
caution is in order. Unavoidably, an aircraft design is 
always a compromise between many contrasting re- 
quirements. As a consequence, the optimum aero- 
dynamic configuration indicated by Eq. (56) is to be 
taken cum grano salis. 


(6) Curvilinear Flight in a Horizontal Plane 


The theory developed for motion in a vertical plane 
can be immediately extended to cover the optimum 
flight conditions in a horizontal plane. It is assumed 
that the tangential acceleration is zero, that the thrust 
is tangent to the flight path, that the velocity is con- 
tained in the plane of symmetry of the aircraft, and 
that the resultant aerodynamic force is also contained 
in the plane of symmetry. 

After projecting the equation of the motion on the 
axes of the principal trihedral associated with the flight 
path (tangent, principal normal, binormal) one obtains 


&, = M, B) — D(x, M. L) = 0 | 
&, =L sin ¢ — (W/g) | = 0> (57) 
®,=Lcosg-—- W=0 | 


where ¢ is the angle of bank and r the radius of turn. 
For simplicity, the problem of extremizing V = 
W (M, L, ¢, 7) is considered, subject to Eqs. (57) and 
to the further constraints 


= B — const. = 0| 
= — const. = 0f (58) 


By analogy with the treatment of Section (3), the 
optimizing condition is written as follows: 


tf The solution to the general problem of optimizing a con- 
figuration which depends upon m arbitrary parameters has been 
indicated by the writer in reference 12. 
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(6.1) Particular Types of Turns 


Three special types of turns are now considered: 
(a) turning flight with maximum angle of bank (VW = 
¢); (b) turning flight with maximum angular velocity 
= Ma(r)/r); (c) turning flight with minimum 
radius of curvature (VW =r). 

The optimizing condition is written in the elegant 
form 


[0 log (7/D)/0 log — 


N(O log D/O log L) [1 — (W?/L*)] = 0 (60) 


where JN is a constant having the following values: 


N=0, forVvV=g 
N =1, for = Ma(z)/r 
N=2, forV=r 


For the particular case where thrust, zero-lift drag 
coefficient, and induced drag coefficient are independent 
of the Mach Number, Eq. (60) can be shown to yield 
the simplified relationships already calculated by the 
writer in reference 8. 


(7) Conclusions 


An analytical method is developed for investigating 
extremum problems of the quasi-steady mechanics of 
flight. A generalized solution is obtained in a detri- 
minantal form. This solution unifies all the optimum 
conditions of motion in a vertical plane into a single 
equation. An analogous result is obtained for the 
optimum conditions in a horizontal plane. The appli- 
cation of the above method to design problems is pointed 
out. 
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(8) Remark 


The cases of flight in a vertical plane and in a hori- 
zontal plane are, actually, particular aspects of a more 
general flight condition—namely, helicoidal flight. 
In this connection, the writer has developed another 
general determinantal equation, valid for optimizing a 
helicoidal flight condition. Considerations of space, 
however, have inclined the writer to omit the treatment 
of the helicoidal case from the present article. 
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On the Optimum Design of an |-Section Beam* 


S. Krishnan and K. V. Shetty 
National Physical Laboratory, New Delhi, India 
March 16, 1959 


SYMBOLS 


= external moment 
2h = web depth 
= flange width 
k = ratio h/b 
t = wall thickness 
A = area of material in the section 
I = moment of inertia 
K = constant for local buckling stress 
oy = yield stress 
E = Young’s modulus 
Lagrange multiplier 


Db pocagennn DESIGN of any thin-walled beam is generally based 
on the condition of simultaneous failure in local buckling 
of the compression element and by direct stress due to bending. 
A method of determining the optimum dimensions for an I-sec- 
tion beam by minimization of the weight subject to the above 
conditions of failure is described and the results are compared 
with those obtained by Shanley and Micks! by their well-known 
structural-index method for a typical example. 

Fig. 1 represents an I-section with dimensions as shown. While 
computing the moment of inertia only central-line dimensions 
are used, which is justified for a thin-walled section. It is also 
assumed that web and flange thicknesses are the same. 

The two conditions governing the beam failure may be ex- 
pressed as: 


* The work described in this note is a part of the research program of the 
National Physical Laboratory and is published with the permission of the 
Director. 
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2h 
t 
7 
Fic. 1. Beam cross section. 


KE (t/b)? = oy (1) 
KE (t/b)? = Mh/I (2) 


where J = (2/3)th® + 4bth?. 
for I, Eq. (2) becomes 


Putting h = bk and substituting 


KE (t/b)? = M/[2tb?k(2 + k/3)] (3) 
Area A of cross section is given by 
A = 2bt(2 + k) (4) 


For optimum design, A should be a minimum subject to the con- 
ditions (1) and (3). 

For minimization of the area, ¢ and k are treated as variables, 
and the value of b is obtained by substituting for ¢ in Eq. (1) 
which may be put into the form 

b = [KE/o,]"t (5) 
Eq. (3) can be expressed as 
d(t, k) = 2k (2 + k/3) —-C=0 (6) 
Area A is minimized introducing the restraining condition (6) 
with a Lagrangian multiplier, \. 
6(A + Ad) = O 
Differentiating with respect to ¢ and k, 
(0A /dt) + = 0) (7) 
(0A /dk) + = Of 
Substituting for the derivatives of A and ¢ in Eas. (7), 
2b (2 + k) + d 30k (2 + k/3) = 0) 
2bt + At(2 + 2k/3) = 0 

Eliminating \ from Eqs. (8) and solving for ¢ and k from the 

resulting Eq. together with Eq. (3), 
t = 0.541 [M/KE]}/3\ (9) 
k=1.3 
and substituting for ¢ in Eq. (5) 
b = 0.541 (KE)"6 M189, 

It may be noted that while optimum ¢ and b depend on the 
applied bending moment, & is independent of the moment and is 
the same for any I-section of uniform thickness. 

Applying this analysis to the example worked out by Shanley 
and Micks, optimum dimensions are computed for a maximum 
bending moment of 36,000 Ibs. in. For the conditions one side 
free and three sides simply supported K may be taken as 0.45 


while Z for Al-alloy is 10.8 X 108 psi. Using an yield stress of 
55,000 psi. for 75S-T Al-alloy, values of t, b, and h are com- 


(8) 
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TABLE 1 
Particulars for Optimum I Section 


Flange Web 
Thickness width depth Area Extreme fiber 
in. in. in k in.? stress psi* 
by minimization 0.107 1.97 2.56 1.3 0.495 55,000 
by structural 
index 0.138 1.83 1.85 1.0 0.758 67,000 


* Extreme fiber stresses have been computed on the basis of linear dis 
tribulion. 


puted. The results are shown in Table 1 along with those ob- 
tained by the structural-index method. 

It is seen from the table that for the same conditions of loading 
there can be a saving of about 8 per cent in weight by the method 
of minimization. It may further be noted that whereas the 
minimization method is based on the yield-stress criterion, the 
structural-index method employs the modulus of rupture, which 
is considerably higher than the yield stress. Also the ratio 13 
of web depth to flange width is considered quite reasonable for 
practical design. 
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Steady-State Circular Pitching and Yawing 
Motion of Symmetric Missiles 


Charles H. Murphy, Jr. 
Ballistic Research Laboratories, Aberdeen Proving Ground, Md. 
March 17, 1959 


T IS SHOWN in reference 1 that a spinning symmetric missile 
acted on by certain nonlinear aerodynamic moments will 
eventually perform a steady-state circular motion of its axis 
about its trajectory. It is further shown that such limit circular 
motion cannot be associated with a nonspinning symmetric mis- 
sile whose damping-moment coefficients, Cy;, and Cyq, are linear 
functions of its squared total angle of attack, a? + 8%. (Under 
such conditions, a planar limit motion, however, is possible.) In 
reference 2, it is proved that even if the static-moment coefficient, 
Cva, as well as the damping-moment coefficients are expanded as 
power series in a? + 8?, limit circular motion is impossible for a 
nonspinning symmetric missile. 

Limit circular motions, however, have been experimentally ob- 
served during ballistic range tests of simple statically stable non- 
spinning bodies of revolution (Fig. 1). It is the purpose of this 
note to indicate a possible moment nonlinearity which could 
produce this observed motion. 

If we assume an essentially straight trajectory, the derivatives 
of the angles of attack and side slip may be related to the angular 
velocities—i.e., 


(1) 


Next, it is assumed that the aerodynamic moment can be ex- 
panded as power series in a, 8, &, 8, g, and 7. Using results of 
Maple and Synge,’ it is shown in reference 4 that all the terms in 
this power series, which remain when the requirement of complete 
rotational symmetry is used, can be written in the following form: 


a = 


M + iN = + (Cu, + Curae’) (2) 
Core = EE)" EE)" E'E)? (3) 
Cu, + Cua = EE)™( EE)" (4) 


where — = B + ia, &’ = il/V = (B + ia)l/V, m,n, and p are 
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nonnegative integers and primes denote derivatives with respect 
to dimensionless arclength s = f(V/I)dt. 

In references 1 and 2, » and p were forced to be zero. We will 
now consider all zeroth and first-order terms in Eqs. (3) and (4). 
The difierential equation for — assumes the form 


+ He’ — Mt =0 (5) 


where = (pSl‘2m)(Cya — — + Cura)] 
= How + + How + Hooé’t’ 
M = 
= Mow + + Mowti’ + Moné’é’ 
= mass 
ky = NI ml? = V 1,/ml? 


Under the quasi-linear assumptions,! the actual nonlinear 
angular motion is approximated by a pair of rotating exponen- 
tially damped two-dimensional vectors whose frequencies and 
damping exponents are functions of amplitude. 

Therefore, = Kie'* + (6) 
where $; = + 
K;’ = —Kj; 
and @’;, \; are functions of the K; 
For the case in question, 
2 [Hoo 1? + Mow( Ki? = 2K,") 

MocoHon( Ki? + (7) 
and a similar expression for A». Since Mo) is usually not bigger 
than 10° 4, the last term in Eq. (7) can probably be neglected. 

Amplitude planes can now be constructed by the use of the 
equation 

d( d( = Ko?) (8) 
This equation predicts a limit circular motion with magnitude 
Ho/(Ai00 + Mo) 
when (a) Hoy < 0 
(b) — Mow < < —38 Mao 

The amplitude plane for (Hio0/How) = —200, Mo:o/How = 100 
is shown in Fig. 2. The actual & motion for these values of 
Hywo/ Hoe and Mow/Heoe and Mooo/Hoe = 0.003 was computed 
from Eq. (5) and is that shown in Fig. 1. 

Thus, nonlinear moments of the form and 
will produce limit circular motions. Estimates‘ of these terms 
for a slender body can be made by inserting the effect of angular 
velocity on the cross velocity in expressions for the nonlinear 
Limit cycles have been observed 
In order 


moment given in reference 5. 
for low-fineness-ratio bodies for which 9 is positive. 
to determine whether the moment terms of this note actually 
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Fic. 2. 


do cause the observed limit cycles, fairly difficult direct measure- 
ment would have to be made. 
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On the Effect of Flow Discontinuities on 
Oscillatory-Pressure Measurements in 
Supersonic Flow* 


Erik Mollo-Christensen and John R. Martuccelli 
Associate Professor and Senior Research Engineer, Respectively, 
Massachusetts Institute of Technology, Cambridge, Mass. 


March 18, 1959 


N THE COURSE of pressure-distribution measurements on oscil- 
latory wedges in a supersonic tunnel, it was found that rela- 

tively slight flow discontinuities can cause large errors in the 
measured oscillatory pressures. 

References 1 and 2 give the results of an experimental investi- 
gation of the effects of thickness, mean angle of attack, and 
amplitude of oscillation on the pressure distribution on a wedge 
oscillating in supersonic flow. The measurements were per- 
formed on a wedge of 6.88° included angle at Mach Numbers 
1.4 and 1.8. The pressures were recorded using piezoelectric 
pressure transducers connected with the wing surface through 
small pressure conduits. 

During the course of the investigation, it was discovered that 
the oscillatory pressures on the wing surface varied from one 
month to another under otherwise identical conditions. After 
all other possible sources of error had been eliminated, an investi- 
gation was made to determine the effect of flow discontinuities 
on the measured pressures. 

Consider Fig. 1. A flow discontinuity intersects the leading- 

* This work was in part supported by the U.S. Navy, Bureau of Aero- 
nautics under Contract No. NOa 58-313-d. 
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TUNNEL WALL 


LE SHOCK 


TUNNEL WALL 


Fic. 1. Interaction of flow discontinuity and oscillating wedge. 


edge shock at 4 and the wedge at B. As the wedge oscillates, 
the leading-edge shock will oscillate, and the change in refraction 
conditions at A will make point B move back and forth over the 
wedge surface. If there should happen to be a pressure hole in 
the path of B, the oscillatory pressure would change from that 
which would be observed in the absence of the flow discontinuity. 

To estimate the magnitude of the pressure signal caused by the 
presence of the discontinuity, consider a Mach-wave discon- 
tinuity. One may then assume isentropic conditions across it, 
and one has 


= {1 + [(y — (1) 


where Py = stagnation pressure, P = local pressure, 1 = Mach 
Number, and y = ratio of specific heats. Since P» is constant, 
differentiating Eq. (1) yields 


dP = —(yMP/{1 + [(y — 1)/2]M%})dM (2) 


which gives the pressure change due to a change in M. 

For dM = 1/100, which is near the limit of resolution of a 
schlieren system, one finds under the particular conditions of 
our experiments that dP is of the same order of magnitude as the 
oscillatory pressures which would be present in the absence of 
the flow nonuniformity. The effect of these flow discontinuities 
has been verified by a series of tests in which flow discontinuities 
were intentionally introduced into the flow. 

The obvious way of reducing such errors is to increase the 
amplitude of oscillation. This is, however, difficult due to 
power limitations in the oscillator mechanism. Another solution 
would be to measure the pressure over an area larger than that 
traversed by a flow discontinuity. 

It is felt that this may be the explanation for the lack of pub- 
lished experimental data on pressure distributions on oscillating 
wings in supersonic flow. 
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Analysis of Effects of Diffusion of a Foreign Gas 
Into the Laminar Boundary Layer of a 
Supersonic Flow of Air in a Tube 


John R. Radbill and Joseph Kaye 


Research Engineer, Aerojet-General Corporation, Azusa, Calif., 
and Professor of Mechanical Engineering, Massachusetts Institute 
of Technology, Cambridge, Mass., Respectively 


April 8, 1959 


HE PROBLEM of diffusion of a foreign gas by means of uniform 
injection into the laminar boundary layer of a supersonic 
flow of air, with known and controllable boundary conditions, has 


been studied experimentally and analytically in the past 3 years 
in the Mechanical Engineering Department of Massachusetts 
Institute of Technology. This program has been financially 
supported by the Office of Naval Research and the Air Force 
Office of Scientific Research.* This note gives a brief summary 
of our theoretical studies, taken from reference 1. 

The model chosen corresponds to supersonic flow of air in a 
tube, with a laminar boundary layer present, and with controlled 
uniform injection of a foreign gas, such as helium, for a finite 
length of the entrance flow. The flow model is shown sche- 
matically in Fig. 1. 

The basic equations were derived by means of kinetic theory, 
statistical mechanics, and the thermodynamics of irreversible 
processes, and then were expanded into cylindrical coordinates 
and subjected to an order-of-magnitude analysis to yield bound- 
ary-layer equations. Thermal diffusion and the inverse diffu. 
sion-thermal effect were neglected. The following boundary- 
layer equations apply: 

(0/07) (pur) + (0/02) (pe?) = O (1) 

Momentum: [pa(dw/07) + pw(dw/0Z)] = 
(2/Repn (0/07) [a7(OW/OF)| — (1/kMo?) (Op/0F) (2) 
op/or = 0 (3) 

Single species continuity: (ScoRep,/2) [pa(dc,/07) + 

= (1/F) (0/2F) (4) 

Energy: (Rep,/2)pZp,[Ci(Cpyp — 1) + 1) X 
+ wW(OT/dz)] = (1/Pre) (1/7%)(0/dF) + 
(ADi2/ Seo) (0e1/07) (0/07) { Tep( — 1)} + 
[Rep(k — 1)/2k] W(dp/dz) + My%(k — (5) 


Continuity: 


State: p = pRT (6) 
The boundary conditions are: 
=0, pi = OT/27 = 0 
r7=1 (7a), (7b), (7c), (7d) 


= (Rep /2)Scol — 1) (pu 

7=0: T=T. CG =0 (8a), (8b), (8c) 

where (7), is a constant, and #, and 7, are determined by a one- 

dimensional analysis of the isentropic core. 

A stream function defined by 

pur = —(O0W/03), (9a), (9b) 

is introduced to satisfy Eq. (1) identically, and a coordinate trans- 

formation defined by 


= (1 — 


pir = ow/d7 


(10a), (10b) 
¢ = 2[(L/D) (1/Rep,)]'!2 = (2%/Rep,)'/? 


is executed. Convergent power series in ¢ with coefficients which 
are functions of 7 are assumed to exist for dependent variables in 
the boundary layer, except for pressure which is a function of ¢ 
only. 

The resulting ordinary differential equations were integrated on 
the IBM 704 digital computer at the M.I.T. Computation 
Center. 

The solution, shown in Fig. 2, of the equations for uniform 
injection of helium into a laminar supersonic boundary layer 


* Reproduction of this material is permitted for any purpose by the 
United States Government. 
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along tube length for various values of the injection Reynolds 
Number, Rep; = —Rep (pt)wan/(pW)o. 


with an adverse pressure gradient in a tube, with an entrance 
Mach Number of 5, predicts a small rise in recovery factor, as 
shown in Fig. 2. At tube entrance, corresponding to the lead- 
ing edge of a flat plate, the predicted increase in recovery factor 
is 0; slightly downstream, this increase varies from 0 to about 8 
per cenit at a value of ¢ of 0.035 where validity of the solution 
becomes questionable. 

A parallel experimental investigation of diffusion of helium 
into the laminar boundary layer of a supersonic flow of air in the 
entrance region of a tube has been under way for 3 years. Some 
experimental results? for uniform injection of helium into an 
air stream with an inlet Mach Number of about 4.8 and an inlet 
diameter Reynolds Number of 123,000 confirm the theoretical 
predictions given above. These experimental data show prac- 
tically zero increase in recovery factor ratio with increasing values 
of injection Reynolds Number Rep,;, as long as the flow in the 
boundary layer remains laminar during injection. 
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A Semiempirical Relation for Laminar 
Separation* 


James L. Amick 

Associate Research Engineer, The University of Michigan, 
Ann Arbor, Mich. 

March 23, 1959 


A THEORETICAL RELATION for the pressure rise to a laminar 

separation point in two-dimensional supersonic flow has 
been derived by Gadd.? For separation from a flat surface, 
this relation can be expressed as 

(p, — ~ — 1)Rei)'”* 

where Re; = pi Vix1/ui, x1 is the distance from leading edge to 
beginning of boundary-layer interaction, and the subscript 1 
refers to conditions just ahead of the interaction and outside 
the boundary layer. This equation can be simplified by sub- 
stituting 1,2 for M2 — 1. (For M, = 1.5, this simplification 
changes the pressure rise by 16 per cent; for M, = 2.5, by only 
4 per cent.) There results 


* This work was supported by the U.S. Air Force as part of Contract No. 
AF 33(616)-3573, monitored by Aeronautical Research Laboratory, WCLJD, 


WADC. A complete account is given in reference 1. 


READERS’ 


603 


FORUM 


— ~ Vx 
where x: is the hypersonic viscous-interaction parameter 
x = Re, 


(This parameter has been found to be important in the problem 
of the pressure rise caused by turning of a hypersonic flow to 


allow for thickening of alaminar boundary layer. Its appearance 
in the laminar separation equation is not unreasonable since 
this also is a problem involving turning of a supersonic flow to 
follow a viscous layer.) 

Experimental data on separated flows with completely laminar 
boundary and mixing layers have been obtained by Chapman, 
Kuehn, and Larson.* The values they obtained for the plateau 
pressure rise after laminar separation from a flat surface (plotted 
in Fig. 1 as a function of the hypersonic interaction parameter, 
xi) are well represented by the semiempirical equation 


— = 127V (1) 


for Mach Numbers above 1.3. It is seen that the plateau 


pressure rise varies linearly with V x, in agreement with the 
approximate equation for pressure rise to the separation point 
derived above from the theory of Gadd. 

Eq. (1) is essentially a local relationship between the plateau 
pressure rise and the location of separation for given free-stream 
conditions. In order to solve many problems of pure laminar 
separation, a second relationship between these two variables 
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must be found. In the simple case of separation ahead of a step, 
the reattachment point is fixed, and the oblique-shock theory 
gives the required second relation between separation-point 
location and pressure rise. For most other problems, however, 
the reattachment point is not known. It is, therefore, necessary 
to find a relationship between reattachment-point location and 
pressure rise at reattachment. Further experimental work is 
needed to define such a relation. 

An expression similar to the present semiempirical relation 
has been proposed by Guman.‘ His expression, however, uses 
a Reynolds Number based on the distance to the point where 
the disturbance would intersect the boundary layer if there 
were no interaction. It, therefore, gives the pressure rise for 
given flow conditions directly without requiring local expressions 
for separation-point or reattachment-point locations. An 
empirical constant is used to fit the equation to experimental 
data for a given body shape. Guman’s equation can be written 


— = — x, (2) 


where C is the empirical constant. A comparison of Eqs. (1) 
and (2) is shown in Fig. 2 for the case of a forward-facing step 
at M, = 3.0. 
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Heat Transfer to a Vaporizing, Ablating Surface 


Richard G. Fleddermann* 
Senior Engineering Scientist, Radio Corporation of America, 
Missile Electronics and Controls Department, Burlington, Mass. 


March 23, 1959 


HE DETERMINATION of the heat transfer to an ablating 

surface is an extremely complex one, which has been treated 
theoretically by many authors—e.g., Lees.!. The exact solution 
for an arbitrary material involves knowledge, available for very 
few materials, of surface temperature, surface heat of reaction, 
vapor species, reaction rates of species with boundary-layer 
gas constituents, transport properties of vapor, gas constituents 
at the edge of the boundary layer, and many other properties. 
However, the results of available theories for injection of simple 
gases point the way to a simpie engineering approach as well as 
a straightforward experiment for determining the quantities 
needed for this approach. 

Implicitly bound up in reference 1 is a relation between gq, 
the energy-transfer rate per unit area to a surface with mass 
addition, and qo, the zero mass addition heat-transfer rate, 
where both q and qc are determined for the same surface tempera- 
ture, 7,. This relation is 


q/qo = f(B) (1) 
where 
B my/peteC 
My = mass rate of vapor injection per unit area 


* Formerly, Chairman, Research Branch, Aerodynamics Section, Research 
and Advanced Development Division, AVCO Manufacturing Company. 
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Pe, Ue, He = gas density, velocity and specific stagnation 
enthalpy at edge of boundary layer, respec. 
tively 

gas specific enthalpy at wall with no mass jn. 
jection 


Huw 


This is the same general form used by Rubesin and Pappas? 
for gas injection into a turbulent boundary layer on a flat plate. 
They also used a similar formula for shear foree—namely, 


t/ro = fil Bi) (2) 


where = m2/(peteC 


It is obvious from physical considerations that f(0) = fi(O) = 1, 
Eq. (1), after substitution of Ci, can be expanded to give 


= 1 — n[m.(He — Hwo)]/go +... + (3) 
or neglecting higher-order terms 
qd = Qo — HH. — Huo)m, (4) 


Thus, the term »(H. — Hwo)my gives the lowering of energy- 
transfer rate to the wall due to mass injection rate, m,. A 
similar relation for shear stress is 


T = T) — NtleMy (5) 


Scala and Sutton,* made calculations for air injection into air 
with Sc = Pr = 1 over a wide range of wall enthalpies. The 
author found‘ that these calculations could be correlated by 
7 = 0.8, ne = 0.67 over a very wide range of injection rates, 
indicating that the neglecting of higher-order terms was essen- 
tially correct. Baron,’ made extensive calculations at various 
Schmidt and Prandtl Numbers for light gases as well as air. 
Bade,® found that these results could be put into the form of 
Eq. (4) and thatn ~ M-'/? where M is the molecular weight of the 
injected vapor. The variation is probably a function of specific 
heat as well as molecular weight for complex molecules. For 
this data air = 0.74. The available calculations for injection 
into a turbulent flow indicate that 7 is of the order of 0.4 to 0.5 
of the laminar case. 

The energy transfer rate into the wall is given by 


qu = 9 — Gr — (6) 
where g, = net rate at which energy is radiated away from the 
surface 


A» latent heat of vaporization to produce vapor species 


actually present at surface 


The above holds for either solid or liquid surfaces. 

Most of the complexities of the boundary-layer flow can be 
lumped into the one constant n. The simple rules set forth for 
estimation of 7 are difficult to apply for arbitrary wall surfaces 
due to uncertainty in knowing the vapor species at the surface. 
However, the » for laminar flow may be found experimentally 
by considering steady-state ablation at a blunt-body stagnation 
point. At extended times under constant external conditions 
an ablating surface reaches a condition where the mass rate of 
ablation, heat-transfer rates, and temperature distribution in 
the material are constant with time. This is the condition of 
steady-state ablation. 

For this case qu is 


Tw 
qw = mf = — T°) (7) 


for a subliming solid, 


Tw 
= me cdT = — T°) (8) 
T 
for a vaporizing glassy material, and 


Tw TM 
de = [ef + XAxy cat | (9) 
TM T° 


for a solid with a discrete melting point, Ty. T° represents the 
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initial solid temperature and c the specific heat of the material. 
The mass rate of ablation, m, must be distinguished from the 
mass rate of vaporization, m,, for materials that have a liquid 
phase; m™ = m, only for a subliming solid. 

Materials with discrete melting points will be neglected in the 
discussion that follows, as they are usually characterized by low 
coefficients of viscosity, so that the melt is swept away without 
reaching vaporization temperature (thus, m, = 0). The con- 
stant E can be loosely termed the heat-capacity efficiency of 
the liquid layer. It represents the fact that not all of the melt 
reaches the wall temperature, 7,. By specializing to the 
stagnation point, one is assured that 0 Z E Z 1 as some liquid 
will flow away from this point. The above limits may not be 
true elsewhere. 

Combining Eqs. (4), (6), and (7), one gets an energy-transfer 
balance for a subliming solid (m = m,), 


= m(n(H, wo) +A» + tA (10) 
A similar equation for a vaporizing, glassy material is 
go — Gr = m,[n( He — Hw) + rv] + — T°) (11) 


The effective heat of ablation is defined by g* = (qo — gr)/m. 
This gives 


q* = n(H, = Huw) + Av + (12) 


and 


q* = [n( He — Hwo) + rv] (me/m) + Eca(T. — T°) (13) 


for subliming and glassy materials, respectively. 

If a sample of a subliming material is placed in an extremely 
hot air source, such as a plasma jet, the sample willablate and, 
eventually, reach steady-state ablation. The recession rate of 
the surface at the stagnation point can be measured photo- 
graphically and m found by multiplying this rate by the material 
density. gq, can be found by a specialized pyrometer and qo 
by a calorimeter of the same shape as the model. Actually, the 
calorimeter reading must be corrected for an estimated wall 
temperature. Thus, g* can be calculated, as qo, g,, and m are 

In Eg. (12), the wall temperature appears as caTy — 
This term is small compared with A, and 7H,. Thus, 
By run- 


known. 
q* is independent of 7, but is a linear function of H,. 
ning a set of samples over a range of H, values, this linear relation 
can be found. The slope of the curve is 7 and the intercept is 
+ Tw — T°) — nHwo. As Ca can be obtained independently 
of this experiment and 7, estimated from radiation studies, this 
experiment yields both y and X, as well as q,. 

The exact variation of m,/m and E for a liquid cannot be 
determined without extensive calculation. However, in reference 
7 it is shown that they approach unity at very high values of 
H,. Thus, extrapolation of the g*, H, data at high H,, gives 
» and A». This analysis could also hold for a system such as 
graphite-air that has extreme energy release by reaction in the 
boundary layer, provided it is physically possible to obtain 
H, > Hw in the hot-jet facility used in the tests, as the q* 
telation can be shown to be 


q* = — Heo) + 
[ky + Tw — T?)I/{1 + [Qc/(He — Hwo)}} 

A simple experimental solution for turbulent flow is not obvious 
at the present time. For computational purposes, one can use 
the laminar A, and the approximation that the turbulent 7 is 0.4 
to 0.5 of the laminar 7. 

Thus, for any chosen 7, and m, one can calculate the heat- 
transfer rate to an ablating wall, utilizing Eqs. (4) and (6) once 
the applicable 7 and A, values are known. These constants can 
be estimated from theory for simple vaporizing materials or 
determined experimentally for more complicated substances 
by the method outlined in this paper. 


(14) 
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Forces and Moments on Oscillating Wing-Body 
Combinations at Supersonic Speed* 


Haim Kennet and Holt Ashley 
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a slender midwing-body combination oscillating in a 
supersonic airstream (Fig. 1). The linearized equation for 
the amplitude of the perturbation potential is well known to be 


(1 — + + — + = 0 (1) 


All coordinates and lengths have been made dimensionless by 
division with reference length b (wing root chord, or body length 
in the case of body alone). Velocities have been made dimen- 
sionless with respect to flight speed U, and the perturbation ve- 
locity potential @ = ge” withthe product Ub. The flight Mach 
Number is M and k = wb/U is the reduced frequency of simple 
harmonic motion. The linearized boundary condition on the 
surface of the slender configuration is 

¢, = Nn, + tkn (2) 
gy, being the velocity component in the direction of the outward 
surface normal n. 

A combined process of expanding ¢ in increasing powers of the 
aspect ration parameter o (dimensionless maximum semispan of 
the wing) and the reduced frequency was utilized to determine 
lift and moment due to small translational and pitching oscilla- 
The development follows the original scheme of Adams 


2 


tions. 
and Sears! and its application in the transonic case by Landahl. 


* The investigation was supported by the Office of Scientific Research, 
USAF, under Contract AF18(600)-961. 


Fic. 1. Plan and rear views 
of a typical combination of a 
slender plane wing centered in a 
body of revolution, flying at 
supersonic speed. 
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series, it is found that 


Lex, dy + (D4/Dx*) } 


¥, 2) = 


where 


The first term on the right of Eq. (3) is easily recognized as 
that associated with slender-body theory, and it satisfies Laplace’s 
equation in the crossplane. It constitutes the first term in an 
expansion of the form 


= + + — 1)? In — 1) + 


O[k?M‘4o4 In — 1) (5) 
For a detailed derivation, the reader is referred to references 3 
and 4. There, the complete calculations of ¢“ and ¢ are ob- 
tained by means of conformal transformation from the wing- 
body combination to a slit on the y-z plane. 

In applying the foregoing result, attention has been focused* 
mainly on slow oscillations in pitch and vertical translation of the 
configuration shown here, although other examples, such as 
chordwise bending, higher frequencies of oscillation, and short 
afterbodies extending beyond the trailing edge of the wing, have 
also been treated.‘ 

Analytical expressions will be found in references 3 and 4 for 
stability derivatives such as Cag, Cu,, and of plane 
wings, bodies of revolution, and wing-bodies. Only a fewresults 
of particular interest will be reproduced here. 

The influence of plan-form shape on fixed-axis pitch damping 
Cu, + Cmg which cannot be distinguished using conventional 
slender-body theory, is shown in Fig. 2. The curves indicate the 
advantage of leading-edge convexity and slight disadvantage of 
concavity in the lower supersonic flight regime. The situation 
is reversed above approximately M = 1.7 for the example chosen. 
For the case of the wing alone, the values of Cu,and Cug, the 
stability derivatives which contribute to the damping in pitch, 


| 

14 16 18 20 
MACH NUMBER 
Fic. 2. Theoretical fixed-axis damping in pitch at supersonic 
speed for slender midwing-body combinations, with various wing 
plan forms, wing aspect ratio A = 1.5 in all cases, and various 
diameter-wingspan ratio R, Axis of pitch located at a = 0.6. 
Dashed lines indicate supersonic leading edge. 
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D?/Dx? = (M? — 1)(02/dx*) + 21kM%(0/ox) — k?M?, 
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By inverting the Fourier transform of the doublet distribution suitable for representing this motion and expanding in a power 


1, 0*)dn + (2/2) (7/2) /M? — 1] x 
O*)dy — (2/2) (a/ax) In (x — &) 7, | dt — (</4x) /M 4. _ 2]. 


— Ag, 0* dn | dt + O[(M? — In o — 1) + In — 1) (3) 


F=V(y — 10)? +2 (4) 
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Fic. 3. Comparison between theoretical and experimental 
values of static stability and fixed-axis damping in pitch for the 
configuration shown. (For description of model see reference 6.) 


as given by the present theory, are in complete agreement with 
the corresponding values given by Miles,> when terms of the 
same order of magnitude are retained. 

The comparison between an adaptation of the present theory 
and experimental results obtained by the Canadian National 
Aeronautical Establishmentf is shown in Fig. 3. Other reason- 
ably satisfactory comparisons with measured static and dynamic 
stability derivatives will be found in reference 4. 
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On a Cross-Elasticity Phenomenon in 
Symmetrically Nonhomogeneous Plates 
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HE FOLLOWING is concerned with the stress-strain relations 

for an aeolotropic nonhomogeneous plate that has a middle 
plane of elastic symmetry. It will be shown here that the theory 
of transverse bending of such a plate involves a cross-elasticity 
phenomenon which is important for the development of a plate 
theory that includes transverse shear and normal-stress deforma- 
tions, as will be shown in a forthcoming report. 

Let x, y be the coordinates in the undeflected middle surface 
of the plate and z the thickness coordinate. Denoting by ci; the 
elastic moduli of the plate material one can write the following 
stress-strain relations: 


Or = Cues + Ci2€y + | 
oy = + Corey + Coeyzy (1) 
Try = + + Cooyzy 


and we specify that 


cii(—2) = cii(z) (i,j = 1, 2, 6) (2) 
Assuming the displacement components to read 
u = y)z, v = B,(x, y)z, w = wo(x, y) (3) 


the corresponding strains in Eq. (1) will take the form 
Yey = (Bry + (4) 


Defining stress-couples M,, /,, Mz, in the usual way and writing 


+h/2 
tii = (5) 


€z = 2; ty = Byy 2; 


where h is the thickness of the plate; we have first 


M, = + Ci2Bysy + + By,z) 
M, + C22Bysy + + 
Muy = €.682,2 + Cx6By,y + + By,r) J 


We solve Eq. (6) for Bz, z, By, y and (Bz, y + By, z) and introduce 
the result into Eqs. (4) and (1). In this way, we obtain 


= + MyAz + MzyAs) ) 


oy = + MyAs + ? (7) 
try = 2(MzA; + MyAs + MzyAs) } 


(6) 


where A, = ci2Cie + Cis 
As = + + 
A; = en Cis + C12 Ca cis Coe 
A, = + + coe Cis 
As = + C22 Coe + Cop C26 (8) 
Ag = + C22 Cog + Ces 
A; = copCi2 + 
As = ceCiz + + 
cisCis C26 Coe + 
and 
Cu = (1/A) — e262), = —(1/A) — 
C22 = (1/A) (Gules — Cs = (1/4) (Gi2@26 — CisC22) (9) 
Coe = (1/A) — Cx = —(1/A) — 
Cn Cr | 
A= | Cz | (10) 
1 C26 Coe | 


In a plate theory which includes transverse shear and normal- 
stress deformation, the cross-elasticity terms will be shown 
to be responsible for an increase in the order of plate differential 

* The author wishes to express his appreciation to Professors E. Reissner, 
C. H. Norris, and A. G. H. Dietz of M.I.T., for helpful discussions of the 
subject of this paper. 
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equations and the number of boundary conditions corresponding 
to it, as compared to the order of the system and its number of 
boundary conditions for the homogeneous plate. 

In contrast to this, no such change is encountered in a plate 
theory without consideration of shear and normal-stress deforma- 
tion. 


Remarks on ‘‘On the Structure of Jets From 
Highly Underexpanded Nozzles Into 
Still Air’ 


S. L. Bragg 
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March 30, 1959 


T° A RECENT PAPER, Adamson and Nicholls! described a method 
of calculating the boundary of an underexpanded supersonic 
jet. An even simpler analysis, based in two-dimensional flow, 
and described in detail below, may be applied when only the 
maximum diameter of the boundary and not the wave length is 
required. 

Consider a stationary jet exhausting into surroundings at a 
pressure, Po, lower than the static pressure, Py, in the exit plane 
of the nozzle. Let A7 be the throat area of the nozzle, Ay the 
flow area at the exit, and A,, the maximum cross section of the 
expanded jet. 

We assume that flow conditions are uniform across each of the 
latter two sections. If the corresponding mean flow velocities 
are Vy and V,,, and P,, is the mean pressure in the plane of maxi- 
mum area, the equation of motion of the*mass flow, M, reads 
(MVw/g) + AwPn + (Am — An)Po =(MVn/g)+ AmPm (1- 

Since only small-amplitude waves occur between sections N 
and m, we may consider the flow to be isentropic. By applying 
the equations of constancy of mass and energy and the perfect- 
gas relations, it is easily shown that, for isentropic flow through 
a varying area, A, the momentum parameter, or thrust coefficient 
in vacuo, (MV/g + AP)/ArP; is a function only of the area 
ratio A/Ar-P, is the (constant) stagnation pressure of the flow. 


That is, ((MV/g) + AP]/ArP, = f(A/Ar) (2) 
Eq. (1) then reduces to the form 
f(An/Ar) + ((Am/Ar) — (An/Ar)|Po/Pt = f(Am/Ar) (8) 
This equation for A,,/Ar can be solved graphically on the plot 
of thrust coefficient vs. A/A7, by drawing a straight line of slope 
P./P, through the point representing conditions at N. The 
point where this line cuts the curve again then gives conditions 


for point m. 
Fig. 1 shows how the maximum area of the expanded jet is 


found in this way. 
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of Blasts as Traveling Gusts 
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gestae OF THE COMPLEXITY of predicting the forces on an 
airfoil resulting from penetration of a traveling blast wave, 
a simple approximation that is frequently made is to represent 
the blast wave as a traveling gust, neglecting the associated 
overpressure. The resulting lift-growth function can be con- 
ceived of as resulting from two causes: growth of circulation 
and a noncirculatory part which, in incompressible flow at least, 
can be attributed to apparent mass. The circulatory lift 
reaches a steady-state value whereas the “appare:t mass” 
effect has a finite impulse. Those who have used the traveling 
gust representation of the blast wave have invariably been aware 
of the associated overpressure and the diffraction of the blast 
wave about the airfoil. The force associated with the undif- 
fracted overpressure is proportional to the area of the airfoil 
subtended by the blast at any instant, and thus yields a finite 
impulse which is proportional to the volume. For a thin airfoil 
this impulse is negligible. The force associated with the dif- 
fraction of the blast wave, on the other hand, is not negligible, 
and it is not clear how much of this diffraction is accounted for 
by the traveling-gust representation. 

It is the purpose of this note to show that when a blast wave of 
any magnitude is represented by a traveling gust and the im- 
pulse is computed using apparent-mass concepts, the result is 
equivalent to accounting for diffraction of the blast wave. 
Since the familiar lift-growth functions for a traveling gust in- 
clude apparent-mass effects as well as the circulatory lift growth, 
these functions are applicable directly to the blast-penetration 
problem because the diffractive loading has already been taken 
into account. 

Criscione and Hobbs,' using a simplified model based on shock- 
tube investigations estimate the diffractive loading on a thin 
airfoil by a shock wave (see Fig. 1). Inasmuch as the airfoil is 
thin and stationary, forces both due to undiffracted overpressure 
and growth of circulation are absent. The model of Criscione 
and Hobbs has subsequently been partially justified theoretically 
by Ehlers and Shoemaker.? Assuming for simplicity that all 
shock and rarefaction waves travel at the speed of the incident 
wave (an assumption slightly different from that made in refer- 
ence 1), the average pressure differential is 


p(t) = po)/2] + fi {((8 + 6po,)] 
(bo, + 1)/2]} [1 — (U.t/c)] (1) 
where p) = pressure of the ambient air; p; = pressure behind the 
incident wave; ¢ = time from instant of impingement; ¢ = chord 
of airfoil; U’, = shock velocity, and po, = po/p1. 
Integrating the resulting average force over the time for which 
p(t) is positive, the impulse due to diffraction is 


Ts = (4/3) + 4)?/[(6po, + 1) (6pe, + 5)]} X 
(c?Ap/U;) (2) 
where Ap = pi — po (= overpressure). 
On the other hand, the impulse due to apparent mass of a 
traveling gust striking a thin airfoil ia incompressible flow is 
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Fic. 1. Diffraction of blast by thin airfoil (schematic). 


Iq = (7/4)pw c? (3) 
where p = density of the medium, and, w = material velocity of 
gust. 

During the period of lift development, the airfoil is always in 
the region behind the shock. Hence, it is natural to choose p in 
Eq. (3) as the density behind the shock, p;. In order to com- 
pare Eqs. (2) and (3), it is necessary to insert the following Ran- 
kine-Hugoniot relation into Eq. (2): 


Ap = = pil(6pq + 1)/(po, + 6)] (y = 1.4) 
Eq. (2) becomes 
Ta = f(po,)pi (4) 
where f(po,) = (4/3) {(3po + + 5) (po, + 6)]} 


This is the same form as Eq. (3) with f(o,) replacing 7/4. 
Examination of f(o,) reveals that f(o,) varies, approximately 
linearly, from 0.848 for an acoustic wave to 0.712 for an infinitely 
strong blast wave. In comparison with 7/4, the maximum 
error is 10 per cent. This excellent agreement shows that the 
impulse on a thin airfoil due to diffraction of a blast wave can 
be accounted for by using the traveling-gust model and com- 
puting the impulse due to apparent mass in a medium of density 
p:. For a moving thin airfoil, Hobbs’ shows that the apparent- 
mass impulse is independent of the forward speed of the airfoil. 
It might then be speculated that the diffractive impulse would, 
likewise, be independent of the forward speed. Therefore, it is 
suggested that functions predicting lift growth based on the 
traveling-gust representation are directly applicable provided, 
merely, that p; is used. 

Whereas only impulses were compared and shown to be identi- 
cal, their force-time behavior may be quite different. _Neverthe- 
less, because the time durations of both the apparent-mass effect 
and the diffractive effect are much less than the time required 
for the circulatory lift to build up, the detailed time history of the 
pulse is usually of no consequence. 
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